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The Elimination of Tuning-Induced 
Burnout and Bias-Circuit 


Oscillations in IMPATT 


Oscillators 


By C. A. BRACKETT 
(Manuscript received August 18, 1972) 


IMPATT diode microwave oscillators suffer from the effects of low- 
frequency instabilities, which include excessive up-conversion of bias- 
circuit noise, bias-circuit oscillations, and diode burnout induced by 
tuning at the microwave frequency. These instabilities are particularly 
troublesome in GaAs diodes, although also present in both Ge and Si to a 
lesser extent. Moreover, these instabilities are more prominent in higher 
efficiency, higher power diodes, presenting a severe systems problem in the 
practical utilization of GaAs diodes at their highest power and efficiency 
levels. In this paper, it is shown that these instabilities may be eliminated 
in a systematic and well controlled manner with litile or no loss tn micro- 
wave power or efficiency. 

It ts shown that the source of the unstable behavior 1s a low-frequency 
RF voltage-induced negative resistance which extends from dc to several 
tens, and perhaps hundreds, of megahertz, depending on the loaded Q of 
the microwave circuit. The negative resistance is an unavoidable fact of 
large-signal avalanche diode operation and 1s due to the rectification prop- 
erties of the nonlinear microwave avalanche. Experiments are performed 
in which the negative resistance and its associated inductive reactance are 
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measured as functions of the de bias current, baseband frequency, and 
microwave circuit loading and Q. 

An analysis is performed and a simple small-signal equivalent circuit 
as derwed for the diode terminal impedance which yields quantitative and 
qualitative understanding of the interaction of the IMPATT device with 
the microwave oscillator circuit. 

Using the equivalent circuit and experimental characterization, a sta- 
bility criterion is developed that is simply applied to any diode-circuit 
combination. Using this criterion, several examples of common configura- 
tions are investigated, including the maximum shunt capacity permissible 
in the bias circuit for stable operation, and the effects of line length in the 
bias circuit. Also discussed are two means of achieving stable operation, 
both of which have been tested experimentally and shown to work. 

The scaling laws are derived in an approximate manner and used to 
show that diodes designed for higher frequencies will have less induced 
negative resistance. Thus, mm-wave silicon IMPATT's would have less 
trouble with low-frequency instabilities than those designed at 6 GHz. 


I. INTRODUCTION 


Low-frequency instabilities including noise, bias-circuit oscillations, 
and diode burnout at low bias currents are often encountered when 
tuning high-power, high-efficiency IMPATT diodes. Gallium arsenide 
IMPATTs, which produce high power at higher efficiency and lower 
noise than either silicon or germanium IMPATTs, are particularly 
prone to these instabilitics. Microwave tuning operations in GaAs 
often burn out diodes at bias currents of one-half to one-quarter their 
thermally expected values. This problem has been so serious that it 
has cast doubt on the practicality of using GaAs IMPATTs at their 
fullest potential power and efficiency. These instabilities also make the 
testing and evaluation of diodes a very tedious procedure and cause 
the operating conditions of the oscillator to be very sensitive to micro- 
wave and bias-circuit load conditions. 

The cause of these problems is an RF voltage-induced negative 
resistance in the IMPATT diode. This negative resistance has a low- 
pass frequency dependence, extending from essentially zero frequency 
up through several tens, and perhaps hundreds, of megahertz. The 
upper frequency limit is determined by the bandwidth of the micro- 
wave circuit. The magnitude and frequency dependence of the induced 
negative resistance are dependent on the microwave circuit tuning. 
It is possible to tune through very high values of negative resistance, 
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even at low de bias currents. If this negative resistance is not stabilized, 
the results can be an excess amount of upconverted microwave noise, 
bias-circuit oscillations, or diode burnout due to the instability achiev- 
ing an excessive amplitude. Stabilization removes the instability and 
eliminates these ‘‘tuning-induced”’ burnouts, as well as the excess noise 
and bias-circuit oscillations. 

In this paper we (2) discuss the physical origins of this negative 
resistance, (77) develop an equivalent circuit which agrees with ex- 
perimental observations, (777) present an experimental characterization 
of the effect, and (iv) discuss the principles and techniques of stabiliza- 
tion which have been used experimentally to remove the instabilities. 

The written history of bias-circuit design is almost nonexistent, 
inasmuch as very little has been known about the low-frequency im- 
pedance of IMPATT diodes. The existence of a low-frequency nega- 
tive resistance effect was observed by Clorfeine and Hughes,! and they 
correctly suggested that this was induced by large RF voltages through 
a rectification effect and postulated that it was the cause of bias-circuit 
oscillation problems. They did not offer any analytic description of the 
effect, nor did they present a systematic characterization of it. 

The construction of bias circuits has been done rather empirically, 
with the understanding that higher-impedance bias circuits give less 
upconverted noise and bias-circuit oscillations. Olson? disclosed the 
use of a resistor in the center conductor of a coaxial cavity to lower 
noise and remove instabilities, and others have used transistor ampli- 
fiers to raise the bias-circuit impedance. Typical examples of bias-circuit 
instabilities are given in an Application Note*® published by Hewlett- 
Packard, along with a discussion of typical bias-circuit designs that 
they have used in silicon IMPATT oscillator circuits. 

The rectification effect itself was contained in Read’s original 
proposal,* and it has appeared in several large-signal analyses since 
then. What has been missing is a simple picture of how the microwave 
circuit completes the feedback loop to provide negative resistance and 
a more definite link between this negative resistance, bias oscillations, 
and tuning-induced burnouts. 

The equivalent circuit and experimental characterization of this 
work permit the elimination of tuning-induced burnout, bias oscil- 
lations, and excessive noise upconversion in a straightforward manner. 

Outlining the paper, in Section II, the link between burnout, bias 
oscillation, and negative resistance is discussed. In Section III, the 
physical origins are discussed and a simple analysis is given, demon- 
strating reasonable agreement with experiment. In Section IV, the 
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experimental characteristics of the total terminal impedance of the 
diode are given as functions of the bias current, baseband frequency, 
and microwave loading conditions. A small-signal equivalent circuit 
of the diode’s terminal impedance is derived in Section V which is 
simple, yet describes in detail the role played by the microwave circuit. 
Stabilization, the elimination of the instabilities, is discussed in Section 
VI, giving the principles, some examples, and possible techniques. 
Approximate scaling rules are derived in Section VII. 

Section VIII presents a summary of the paper and its major con- 
clusions in a form complete enough to give the more casual reader a 
good understanding of the negative resistance mechanism and its 
stabilization. 


II. BURNOUT AND BIAS-CIRCUIT OSCILLATION 


Before proceeding with consideration of the negative resistance 
itself, it may be valuable to establish a link between the existence of a 
negative resistance, bias-circuit oscillations, and tuning-induced burn- 
out. As described later in this paper, a negative resistance has been 
measured and found to exist in the circuit used at frequencies from 
essentially de up to several tens of MHz. Such a negative resistance 
may act as an amplifier of bias-circuit noise, resulting in a modulation 
of the oscillator amplitude and frequency and producing the upcon- 
version of large amounts of noise to the microwave frequency. Under 
certain conditions, the bias circuit may break into sustained oscillations 
producing a combined AM and FM spectrum about the microwave 
signal. When either of the above events occurs (excessive noise up- 
conversion or bias-circuit oscillation), there is grave danger of burning 
out the diode. Specifically, it is found that 


(z) without stabilization of the negative resistance, low-current 

burnouts are frequent, 

(iz) bias-circuit noise and/or oscillations accompany or precede the 
burnout, 

(z2z) stabilization of the negative resistance removes the noise and 
the bias-circuit oscillations, and stops the low-current burnouts, 

(2v) the microwave tuning conditions likely to produce burnouts 
are the same as those that produce the noise and bias 
oscillations. 


Thus, circumstantial evidence strongly links the negative resistance, 
the bias-circuit oscillations, and the low-current tuning-induced burn- 
outs to each other. A specific sequence of events leading to burnout 
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is not postulated. Indeed, many such sequences can be visualized. 
Preliminary scanning electron microscope photographs indicate no 
difference between pure thermal (nonoscillating) burnouts and low- 
current tuning-induced burnouts. Both have the same characteristic 
gold fingers shorting out the junction that were first described by 
Evans, et al.® 


III. PHYSICAL ORIGINS OF THE NEGATIVE RESISTANCE 


The three essential ingredients which cooperate to produce the low- 
frequency negative resistance are 


(¢) a large-signal rectification cffect wherein the dc voltage (at 
constant current) is lowered by an increase in the microwave 
voltage amplitude, 

(ii) the dependence of the diode’s microwave conductance, 
ga(Va,la), upon the microwave voltage amplitude, V., and the 
de bias current, Jz, and 

(227) the microwave circuit constraint 


ga(Vala) + Gw) = 0 


where G(w) is the microwave circuit conductance into which the diode 
oscillates. The first of these, the rectification effect, is summarized by 
Read‘ in the equation (for the Read diode) 


Wr m 
Va ae — TT, 


oA — AWE. Ve + constant (1) 





where 
Va = de voltage, 
W = depletion layer width, 


7 = drift zone transit time, 


e = dielectric constant, 


A = junction area, 


m = (E/a)(da/dE)|»z, = percentage change in the avalanche co- 
efficient a for a percentage change in the 
electric field, and 

E, = critical field for avalanche breakdown. 


Equation (1) is approximate, neglecting higher powers of V.. Studies 
made with the nonlinear program introduced by Blue,® however, in- 
dicate that (1) is an excellent approximation. 
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The net de terminal resistance from (1) is then 


aV ac mV a dV, 
R, = = Es —= (2) 
dla 2WE, dla 
where R,, = Wr/2eA is the space-charge resistance. 
Irom the microwave circuit constraint, we see that the ac voltage 
amplitude V, and the de bias current J, are [for fixed G(w)] not 
independent of each other. In fact, 


Gz, 
aVa Ola Va 


een ey (3) 


a (2) 
OVa Id 


In the usual situation (dga/dIa)v, < 0 and (dga/0Va)1, > 0 (Temem- 
bering that ga is itself negative), so (dV./dIa) > O and the last term 
of (2) represents a negative resistance. It is this negative resistance 
which causes the instabilities in the bias circuit. It is associated, 
basically, with amplitude modulation of the microwave oscillation. The 
rectification property is due to the nonlinearity of the avalanche process 
whereby a sinusoidal field variation about the critical field, L., pro- 
duces many more charges on the positive swing than it does on the 
negative. This means that either the de current increases, or the de 
voltage drops, as the RF voltage increases. 

A simple picture then is as follows. A positive fluctuation of the 
diode current increases the microwave negative conductance (in 
magnitude). This causes the voltage amplitude, V., to increase in 
order to meet the circuit constraint, and the increase in V, requires 
a drop in the de voltage, thereby creating negative resistance. 

Equation (2) can be written 

















m 1 dP out 
2WE.G dla 


where G is the microwave circuit conductance at the oscillation fre- 
quency and Pou is the microwave output power. This form is con- 
venient because G can be estimated and dPout/dla can be easily 
measured in any particular circuit. Doing this, and using published 
data for the other constants, R; has been calculated for typical silicon, 
germanium, and gallium arsenide IMPATT diodes designed to operate 
at 6 GHz. This data is shown in Table IJ and it indicates reasonable 
agreement with the experimental values also shown. Silicon has the 


Ri = Rese 








(4) 
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TABLE I—Puysicau DATA AND NEGATIVE RESISTANCE VALUES FOR 


SILICON, GERMANIUM, AND GALLIUM ARSENIDE IMPATTs 





Si Ge GaAs Source 
C Breakdown (DF) ~0.7 ~0.7 2.5 Measured 
V Breakdown (Volts) 105 45 95 Measured 
W (microns) 5 4 4 Ref. 7, p. 117 
£.(volts/em) 3.5 & 105 2.1 X 105 6.5 X 105 | Ref. 7, p. 117 
E. da = 
San es 5.0 6.1 19.5 Calculated 
dP out 
ars (watts/amp) 13.6 10 24.4 Measured 
G (nhos) 3 X 10-3 3 xX 1073 6 xX 1073 Estimated’ 
R_ (ohms) 65 117 152 
Rs- (ohms) 25 59 25 Calculated 
Theoretical 
R: (ohms) —40 —58 — 127 Calculated 
Experimental 
R; (ohms) —12 —50 —122 Measured 











smallest negative resistance and the largest discrepancy between 
theory and experiment. Germanium has about (—)50 ohms of induced 
negative resistance, indicating that difficulties might be encountered 
if biased through a 50-ohm bias line. The gallium arsenide results 
indicate about (—)120 ohms induced negative resistance. Measured 
values have been everywhere between (—)85 ohms and (—)150 ohms. 
On the basis of this data alone, it would be expected that biasing GaAs 
diodes would be very difficult. 

A complete experimental characterization of the low-frequency 
impedance is given in the next section. 


IV. EXPERIMENTAL CHARACTERIZATION 


The technique used to measure the low-frequency impedance of 
the diode was to induce bias-circuit oscillations and to barely quench 
them by adjusting the bias-circuit impedance. Then the diode was 
removed from the circuit and the impedance of the circuit as seen by 
the diode was measured at the frequency of the bias-circuit oscillation. 
The quenching impedance, Zg, measured in this way, is the negative 
of the small-signal impedance of the diode, Z;, and could be measured 
as a function of the bias-circuit oscillation frequency, the microwave 
loading, and the bias current. The microwave frequency of oscillation 
was always tuned to that for maximum power output. 
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IMPATT DIODE —~ 


Fig. 1—Microwave cavity and bias configuration. The IMPATT diode is mounted 
at the end of a coax line which is coupled to the microwave cavity. The bias circuit 
includes a stabilizing resistor Rs, a resonant circuit Cz, Ls, Re used to induce bias- 
ou Polen and a resistor FR, to isolate the bias circuit from the constant current 
supply Io. 


In order to perform this experiment, it was necessary to use an 
oscillator and bias circuit with the following two properties: (2) a 
bias circuit which could be made stable or unstable at all values of 
current and microwave loading, and (77) microwave loading which 
could be varied continuously from above oscillation threshold down to 
zero external loading. A suitable circuit is one described by Magalhaes 
and Kurokawa?’ and illustrated in Fig. 1. It consists of a waveguide 
microwave cavity coupled to a coaxial transmission line. The wave- 
guide cavity is formed between a movable short and an adjustable iris, 
which is used to couple the cavity to the output waveguide load. The 
adjustable iris is formed by rotating the rectangular waveguide cavity 
with respect to the output rectangular waveguide at a flanged joint. 

The diode is situated at one end of the coaxial line (its distance from 
the waveguide was adjustable for tuning purposes) and the coaxial 
line is terminated on the opposite side of the cavity from the diode in 
a microwave load. The microwave load is an open circuit to de and 
may be represented as a shunt capacitance of about 35 pF at bias- 
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circuit frequencies. The de bias is fed in through the coaxial line 
from a constant-current supply. An isolation resistor of about 200 ohms 
was used to reduce the effects of the parasitic reactances of the power 
supply and its leads on the bias-cireuit impedance. It was found that 
placing a series resistance R, (2-watt carbon resistor) in the bias line 
would stabilize the induced negative resistance formed in the diode. 
A value of R, = 150 ohms was sufficient to stabilize the diode at all 
bias currents up through the thermal limit of the resistor R, and at 
all microwave loading levels. 

The bias oscillation was induced by reducing R, and adding a shunt 
resonant circuit with a quenching resistor Rg. The resonant circuit 


0.7497 1.017 
1.48 — 
9.55~4~ 7.93 2.18 
11.6¢ FREQUENCY, 
015.8 Sade MHz 
©20.6 © 4.55 


@ 28.9 06.94 


©33.2 





Fig. 2—Quenching impedance versus frequency of bias-circuit oscillations for Ge 
and GaAs measured at the diode’s terminals. Bias currents were 100 mA and 120 mA 
for Ge and GaAs respectively. The quenching impedance is the negative of the diode’s 
small-signal impedance, indicating a negative resistance and inductive reactance. 
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was used to tune the bias oscillation to the desired frequency and it 
was brought to the quenched, or small-signal, level by adjusting Ra. 
Measuring the circuit impedance with the diode removed then gave 
the desired measurement of Zg. 

The result of such a procedure at constant bias current and constant 
microwave loading conditions is shown in Fig. 2 for typical 6-GHz 
Ge and GaAs diodes. Since Z, = — Zag, we see from Fig. 2 that the 
diode impedance has a general low-pass characteristic for the negative 
real part, becoming positive at frequencies higher than some frequency 
fm. It is also seen that the reactive part of Z; is inductive. This general 
character to Z; has been observed in all diodes upon which these 
measurements have been made. The frequency fw is the maximum 
frequency at which the bias circuit can be made to oscillate and is 
dependent upon the loaded Q of the microwave circuit, as will be ex- 
plained in Section V. The scatter in the data from smooth contours 
is probably due to the difficulty of obtaining identical microwave 
tuning on successive measurements. The general range of low-frequency 
asymptotes for R, for GaAs was from (—)85 ohms to (—)150 ohms. 
For Ge this range was from (—)35 ohms to (—)75 ohms, although less 
data is available for the Ge diodes. These experiments have been 
attempted on Si 6-GHz IMPATTs as well, with the results indicating 


NET NEGATIVE RESISTANCE Rt ,OHMS 
POWER OUTPUT, WATTS 





50 60 70 80 90 
©,DEGREES — ©,DEGREES — 
MICROWAVE MICROWAVE 
MAX «———— CAVITY ———® MIN MAX <@——— _ CAVITY ————> MIN 
LOADING LOADING 


Fig. 3—Output power and induced net small-signal negative resistance as a func- 
tion of bias current and microwave loading for GaAs. Microwave frequency was 4.7 
GHz, bias-circuit frequency was ~5 MHz. 
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a maximum observed negative resistance of about (—)20 ohms. How- 
ever, with the microwave circuit tuned to best overall operation, these 
6-GHz Si diodes would not sustain bias-circuit oscillation, but only 
emit large amounts of baseband noise as the bias-circuit impedance 
was lowered. This is explained in Section V, where the equivalent 
circuit of this impedance is developed. 

In Fig. 3, the small-signal-induced negative resistance, R,, is shown 
as a function of bias current and microwave loading. This data was 
all taken at a bias-oscillation frequency of ~5 MHz. Also shown is the 
output power. The microwave loading is changed by rotating the out- 
put waveguide relative to the cavity waveguide, the angle between 
them denoted by 0. For 6 — 0, the guides are aligned and the micro- 
wave cavity has its maximum loading, which implies zero (or at least 
minimum) microwave voltage amplitude. Increasing the rotation 0, 
one passes through oscillation threshold to maximum power output 
and finally to maximum oscillation amplitude and zero output power 
in the minimally loaded condition. Thus, 6 may be considered to be an 
indication of microwave voltage amplitude or microwave load con- 
ductance. From Fig. 3 we make the following observations. 


(z) At low bias currents, the negative resistance is small at normal 
maximum power loading conditions, but does become signifi- 
cant at the lightest loading (6 — 90 degrees). 

(iz) The negative resistance increases rapidly with bias current 
until it saturates at a level around (—)120 ohms (for this 
diode). Further increases in current tend only to translate the 
maximum negative resistance toward greater microwave load- 
ing conditions (smaller ©), with the result that, at 160 mA, the 
loadings for maximum power and maximum induced negative 
resistance are nearly coincident. 

(zit) There appears to be considerable structure to the negative 
resistance curves. This structure appears to be real because a 
decrease in either the current or the angle 0 from the condition 
I, = 140 mA, © = 75 degrees results in an increased amplitude 
of oscillation and sometimes burnout of the diode. 


The reason the negative resistance increases with current at low cur- 
rents and increases with 0 at the lower 0 is presumably associated 
with an increase of the voltage V.. The apparent saturation of the 
negative resistance at higher currents and larger 0 is not understood. 
It may be associated with a saturation of the RF voltage V., or more 
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IMPEDANCE 
COORDINATES 


_SMALL—SIGNAL 
IMPEDANCE 





! 
LARGE—SIGNAL 
IMPEDANCE 


Fig. 4—Schematic illustration of the negative of the small-signal and large-signal 
diode impedance at bias-circuit frequencies. The large-signal impedances move in the 
direction of decreasing negative resistance with increasing bias oscillation amplitude. 


likely with the detailed dependence of the diode’s microwave admit- 
tance, ya, on V, and Jq at high V, and high J. 

One further experimental characteristic is that, if the bias oscillation 
is allowed to rise above the quenching level, the real part of Zq@ de- 
creases as shown schematically in Fig. 4. This implies that the device 
impedance Z;, is open-circuit stable and is properly characterized as a 
negative resistance. It also explains the well-known fact that high- 
impedance termination of the bias circuit lowers the AM noise signifi- 
cantly.!° Since the negative resistance acts as an amplifier of noise in 
the bias circuit, the lower the circuit impedance, the higher the gain 
and the larger the noise modulation of the microwave power become. 
In the following section, the small-signal equivalent circuit is worked 
out in detail and most of the observed features are explained on 
theoretical grounds. 

The open-circuit stability and large-signal characteristics have a pro- 
found significance as to the elimination of both bias-circuit oscillations 
and diode burnout, as will be discussed at length in Section VI. 


V. LOW-FREQUENCY EQUIVALENT CIRCUIT 


It is evident by now that the microwave oscillator circuit plays an 
important role in the existence of the low-frequency negative resis- 
tance. The question arises, however, as to the role it plays in determin- 
ing the bandwidth and magnitude of the negative resistance and why 
the reactance associated with R, is always inductive. In this section 
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0 


Fig. 5—Schematic admittance-plane description of the microwave large-signal 
device admittance, ya(Va,Ja), and the negative of the microwave circuit admittance, 
—Y¥(). The intersection of yz with — Y determines the amplitude of the microwave 
oscillation A. and the frequency w. An increase in Ja requires an increase in V, 
from A, to Ag as shown. 


an equivalent circuit is developed, on an analytical basis, which fits 
the data of the previous section, and is used further in Section VI to 
consider the problems of stabilization. 

To develop the equivalent circuit, a quasi-static model is used in 
which it is assumed that the diode’s microwave admittance, ya(Va,la), 
is an instantaneous function of V, and Jq. For example, a sudden change 
in the de bias current is assumed to shift the entire large-signal diode 
characteristic as shown in Fig. 5 from ya(Va,la) to ya(Va, La + Ala). 
Because of the large amount of energy stored in the microwave circuit, 
however, the RF voltage V. cannot change instantaneously. There- 
fore, immediately after the change in bias current, the circuit con- 
straint ya(Va,la) + Y(w) = 0 is not satisfied and can only be satisfied 
by the growth of V, from V, = A, to a new value A; as shown in Fig. 
5. Such growth takes time and the amount of delay is related to the 
bandwidth of the microwave cavity. We can easily imagine that, if Ia 
were to fluctuate rapidly enough, faster than the circuit can respond, 
then V, would maintain its average value and the induced negative 
resistance would disappear. 

Using the quasi-static approach mentioned, the low-frequency 
(baseband) impedance, Z;, of an IMPATT diode is derived in the 
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Appendix. It is found that 














| coe 
Z, = Rs. — ———— (5) 
Ley 
with the definitions 
m sinx dv, 6) 
~ AWE, sin 6| dl, 
and 
y = ———G sin 6. (7) 
OY 
Ow 








Here, ry and s are saturation parameters associated with the large- 
signal behavior of the diode admittance ya as defined in eq. (17) of 
the Appendix, and Y = G+ jB is the admittance of the microwave 
circuit. The angles x and @ are defined in Fig. 6. In the admittance plane 
plot of this figure, — Y (w) is the locus of the negative of the microwave 
circuit admittance in the neighborhood of the oscillation frequency w,; 
the arrow indicates the direction of increasing frequency. The large- 
signal admittance of the diode at frequency w, is assumed to intersect 
the —Y(w) locus at the point P, which in fact determines the oscil- 
lation amplitude, V., and the frequency, w.. The line ya(V.) indicates 





Fig. 6—Illustrating the definitions of the angles 6 and x at the intersection of the 
device admittance ya(Va) and the circuit admittance — Y (w). The intersection angle 
is 6 for microwave voltage variations and is x for bias-current variations. The angle 
6 + x is determined by the device alone. 
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the direction of the admittance plane in which the diode’s admittance 
moves for increased voltage, V., at constant current. The line ya(Jq) 
is likewise the direction in which the diode’s admittance moves for 
increased J, at constant oscillation amplitude, V.. The angles @ and x 
give the intersection angles of yz(Va) and ya(Ia) with the locus — Y (w) 
as shown in the figure. It is well known" that one of the conditions for 
oscillation at the point P is that 0 S 6 S 7, with 6 = 7/2 generally 
giving the best system performance (highest oscillator stability, lowest 
noise). 

From (5) and (6) we see that at very low frequencies (w/y <« 1) 
and if sin xy/sin @ = 1 (a reasonable condition, as will be seen) then 
Z, = R,. — R—- with R_ the same as found on a dc basis in eq. (2). We 
see that the space-charge resistance is a stabilizing factor, tending to 
lower the net negative resistance. At low enough frequencies, thermal 
effects become important and introduce further positive resistance. 
The latter is neglected here due to the relatively low cutoff frequency 
for thermal effects compared to y. 

Considering eq. (6) for R_, we see that, if @ approaches either 0 or 
a, R_ can become very large. In fact, 6 = 7/2 gives a minimum R_, 
everything else being constant Intuitively, one can see why 6 = 0,7 
are critical directions. For these conditions, yz(V.) is nearly parallel 
to —Y(w) and slight changes in the current Jz can cause large changes 
in V, (and the frequency w, as well). 

This dependence on @ is consistent with experimental observations. 
A certain set of Si diodes (6 GHz) were found not to have a large 
enough negative resistance to cause bias-circuit oscillations when 
tuned in what was considered to be the optimum manner. Detuning 
the cavity, however, and readjusting the power output to the same 
value gave very large bias-circuit oscillations with a very definite 
FM microwave spectrum. It is presumed that the FM is mostly due 
to 6 being different from 7/2. 

Separately considered, the angle x could be used to minimize R_ or 
even to make the net resistance positive. The difficulty with this is 
that the total angle 6 + x is determined totally by the diode’s charac- 
teristics and has nothing at all to do with the microwave circuit. The 
large-signal theories of Scharfetter and Gummel and Blue® may be 
used to show that the angle @ + x is close to z over most of the useful 
range of diode parameters (V.,Ja,w.). Comparing Figs. 5 and 6 of the 
paper by Gewartowski and Morris® confirms this conclusion experi- 
mentally, at least for the Si diode they investigated. If 6+ x 7 and 

=~ 7/2, then x 7/2 and sin x/sin 6 ~ 1 for best system require- 
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ments. This dependence of R_ on the angle x indicates, however, that 
in laboratory testing (providing x + @ is not actually 7) it should be 
possible to reduce the value of R_ and inhibit bias-circuit oscillations, 
while at the same time maintaining high power output. Experimentally, 
in coaxial circuits with many degrees of freedom in the tuning it is 
found that burnout and bias oscillations can usually be avoided by 
very carefully increasing the current by small increments and retuning 
slightly. 

The parameter y is shown in the Appendix to be one-half the 3-dB 
modulation bandwidth of the oscillator. y achieves its maximum value 
at 6 = 7/2. We therefore see from eq. (5) that the negative resistance 
has a low-pass characteristic and that Re[Z, — R..] has a cutoff 
frequency w, = y equal to that of the microwave circuit. 

The equivalent circuit of Z, from eq. (5) is shown in Fig. 7. Z; is 
inductive, but the frequency-independent representation involves a 
negative capacitance —C_ = — (1/yR_). We note that, in accordance 
with the cutoff frequency mentioned above, the product (— R_)-(—C_) 
= + 1/y gives the time constant associated with the microwave 
oscillator bandwidth. A Smith chart plot of (—)Z./R,. for several 
values of R_/R,. and with « = w/y as a normalized frequency parame- 
ter is shown in Fig. 8. In this figure we see that the loci of —Z./Rs. 
are straight vertical lines, and that the x = w/y values are obtained 
by drawing straight lines between the points z = 0 — jxand |z| =~. 
For any value of R_/R.-, there exists a maximum frequency for which 
Z, has a negative real part. We call this the maximum frequency of 





R_ 
1+] W/Y 


21=Rgc - 


Fig. 7—Low-frequency equivalent circuit of an oscillating IMPATT diode neglect- 
ing thermal effects. Rs. is the space-charge resistance of the diode, y and R_ are 
defined in the text. 
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Fig. 8—The negative of the equivalent circuit impedance, Z7:/Rs-, normalized to 
the space-charge resistance, versus normalized frequency x = w/y. The straight line 
characteristics of this plot make it useful in estimating the ratio of fa from the 
values of R_, Rsc, and y. 


oscillation, fia, and it is given by 


R_ 
fu = 5" aa (8) 


Thus, the maximum frequency of oscillation may be greater than or 
less than y/27, depending on the ratio R_/R... 

The normalization of Z; to Rs. is convenient for simple construction 
of the loci and finding the ratio 27fx/y, but experimental results are 
more appropriately normalized to the bias line characteristic impedance 
as was done in plotting Zg in the previous section (Fig. 2). There are 
three independent parameters in the equivalent circuit for Z;, namely, 


288 THE BELL SYSTEM TECHNICAL JOURNAL, MARCH 1973 


R_, R;-, and y. By measuring the impedance Z;, at three frequencies, 
one may calculate these parameters. Usually, R,. can be measured in 
dc tests, so that only two measurements of Z;(w) are required. In Fig. 
9, the data for the GaAs diode of Fig. 2 have been replotted along with 
the equivalent circuit curve obtained from the data R.(w = 0) = Re. 
— R_= —110Q, fw217 MHz, and R,, = 33.6 Q. The scatter in 
the data points is, as was mentioned before, due to the difficulty of 
disassembling the microwave circuit and reassembling it to the same 
RF conditions after every measurement. Under the circumstances, the 
agreement is considered to be good. 

To illustrate the meaning of y further, consider the simplest model 





Fig. 9—Comparison of experimental (solid dots) and theoretical (solid curve) 
frequency dependence of the small-signal impedance Z;. GaAs diode, J, = 120 mA, 
fo = 4.7 GHz, output power = 1 watt, R- = 143.60, Rs. = 33.60, and y = 59 X 108 
radians per second. 
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of oscillator and circuit possible in which only the diode’s conductance 
is a function of V, [this puts r = 0 in eq. (7) ] and the circuit conduc- 
tance does not depend upon frequency. Then eq. (7) becomes 


sG SW 6 
Y= = (9) 


(2) ~ 201 

Ow 

where Q, is the loaded cavity Q. Since at maximum output power 
s = 2,8 we can write under conditions for maximum power 








y = 2r-Af (10) 


where Af = the full 3-dB bandwidth of the loaded cavity. Equation 
(8) then becomes 


R_ ithe 
fu = ane 1-2 wae (11) 


This dependence of fw on Qr was tested by making the experimental 
waveguide cavity length ~ 3\,/2 instead of ~,/2 and finding the 
maximum frequency of oscillation. This should raise Qz by a factor 
of three and therefore decrease fy by a factor of one-third. The meas- 
ured result was a decrease in fy by a factor of 0.41. This discrepancy 
may again have been due to not achieving the same RF conditions in 
the two cases. 

In another test of eq. (11), the loaded cavity Q was measured after 
having made the bias-circuit impedance measurements necessary to 
determine fy, R;-, and R_ for the same microwave circuit conditions. 
The results were y = 113.3 X 106, R,. = 385 2, R_ = 119 Q which, 
together with the assumption s = 2, give Q, = 261. The microwave 
measurement yielded Q; & 268, in good agreement. 

We conclude that eq. (5) is a good representation of the terminal 
impedance at baseband frequencies in a large-signal IMPATT diode 
oscillator and that y is well approximated by either (9) or (10). 








VI. BIAS-CIRCUIT STABILIZATION 


Up to this point, the major emphasis has been on characterizing the 
induced negative resistance and discussing its role as the cause of 
bias-circuit oscillations and tuning-induced burnout. We now turn 
our attention toward the principles and techniques of the stabilization 
of this negative resistance. By stabilizing we mean the achievement 
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of small-signal stability so that any fluctuation in bias-circuit current 
or voltage eventually decreases to zero. In addition to small-signal 
stability, we also require large-signal stability so that transients in 
power supply lines cannot cause instabilities to develop. 


6.1 Stability Criterzon 


The stability criterion is derived by starting with the loop impedance 
equation 


(Zp + Z1)tB = €n 


where Zz is the bias-circuit impedance, 7, the fluctuating component 
of current, and en is an assumed noise voltage. In the complex fre- 
quency plane (s = o + jw), if Zp + Z, has any zeros in the right-half 
plane (o¢ > 0), that component will grow in time, and the circuit plus 
diode is unstable. The Nyquist criterion is used to determine the 
number of zeros minus the number of poles of Zz + Z; that are in the 
right-half plane. This is not so convenient here, since we are working 
with Z, and Zp, as two separate sets of data and can only change the 
design of Zz. We therefore consider the function Zg — Zs, where 
Zo = — Z,. This must have the same number of poles and zeros in 
the right-half plane. Now, however, we plot the loci of Zg and Zz 
separately as w goes from — © to + and consider the vector pointing 
from Zp(w) to Za(w). We then see that the net number of rotations 
of the vector Zq(w) — Zz(w) indicates the number of zeros minus the 
number of poles in the right-half plane. In addition, since Z, is a 
passive driving point impedance function, it cannot have any poles or 
zeros in the right-half plane. Therefore, 
ZoQ(s) Z3(s) 1 + a Ris Z3(s) 

has no poles in the right-half plane and the number of rotations of the 
impedance vector Zag(w) — Zz (w) gives directly the number of right- 
half plane zeros. 

To understand this in practical terms, we show schematically in 
Fig. 10 three conditions; stable, unstable, and conditionally stable for 
a typical bias circuit. The bias-circuit impedance, Z(w), and the diode 
quenching impedance, Ze(w), are separately plotted as a function 
of frequency over the complete frequency range for which Z;(w) has 
negative real part. The criterion then becomes: The bias circuit is 
stable if, at the point P where the two loci intersect, the frequency 
wpa on the Zg locus is lower than the frequency wpz on the Zz locus; 
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Fig. 10—TIllustration of (a) stable, (b) unstable, and (c) conditionally stable bias- 
circuit impedance loci. 


otherwise, an instability exists. In Fig. 10a, the criterion is satisfied, 
and the circuit is stable. In Fig. 10b, the criterion is not satisfied, 
and the circuit is not stable. In Fig. 10c, there is not even an inter- 
section of Zz with the small-signal locus Zg. However, Zz(w) does 
lie in the large-signal region of Zg (see Fig. 4), and such a condi- 
tion is conditionally stable. At small-signal levels, no growing root 
appears, but if the excitation becomes large through some disturbance, 
then a large-signal instability can exist. We note that changes in the 
bias current or microwave loading can move Zg all the way down to 
zero impedance and therefore such a circuit configuration as shown in 
Fig. 10c would be almost certain to burn out the diode. 


6.2 Some Examples 


We now give some simple examples which illustrate the above 
criterion and shed light on some commonly used biasing schemes. 
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R_—Rsc Ro 





fu 


(a) (b) 


Fig. 11—Stability considerations for constant-current power supply having a 
lumped shunt capacitance C’,.: Section 6.2.1. 


6.2.1 Constant-Current Supply—Maximum Equivalent Lumped Ca- 
pacitance 

The circuit of Fig. lla represents a constant-current supply with 
large shunt resistance and perhaps large shunt capacitance. Part (b) 
of the figure shows the impedance diagram assuming that R,—- ©. 
The maximum capacitance Co max that can be tolerated before inducing 
instability is that which has a reactance equal to that of Zg at f = fm, 
the maximum oscillation frequency. This is 


1 1 
y R_ ao fixe 


Co max = 


Rs 





(a) (0) 


Fig. 12—Stability considerations for constant-current supply with added series 
resistance, R;: Section 6.2.2. 
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(a) (b) 


Fig. 18—Stability considerations for constant-voltage supply: Section 6.2.3. 


For the GaAs oscillators used for the data in Figs. 2 and 3, we 
had R_ — R,. 120 ohms and vy & 2r X 10’ rad/s, which gives 
Co max = 140 pF. 


6.2.2 Constant-Current Supply—sSeries Resistance 


The addition of series resistance between the shunt capacitance and 
the diode, as in Fig. 12a, gives the impedance diagram shown in Fig. 
12b. Stability is predicted for R, > R_ — Rec. It should be noted here 
that this is the easiest possible way to stabilize the bias circuit, and 
were it not for the large de power dissipation incurred in R,, this would 
be ideal. The larger R, is, the larger the difference between Zz (w) 
and Zg(w), and the less the noise amplification is in the bias circuit, 
giving less upconverted noise in the microwave spectrum. 


6.2.3 Constant-V oltage Supply 


For a constant-voltage supply, the series resistance is usually very 
small, and C> is large. Figure 13 shows the circuit and impedance 
diagram and we see that this is a conditionally stable circuit. As the 
bias current is turned up from zero, the Zg locus moves through the 
low-impedance region, and this scheme is certain to burn out diodes. 
The addition of a large series resistance moves the whole Z» locus to 
the high-impedance region and stability is indicated again. 


6.2.4 Constant-Current Supply—Series Inductance 


The addition of a large series inductance, as shown in Fig. 14, might 
be thought to hold the current more constant and therefore add sta- 
bility. In fact, such an addition creates a series resonance which lowers 
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Ls 





(a) (b) 


Fig. 14—Stability considerations for constant-current supply with added series 
inductance: Section 6.2.4. 


the bias-circuit impedance at all frequencies below the resonance. 
Thus, the Z,(w) locus must intersect the Zg(w) locus at a lower fre- 
quency (on Zz), and this is in the direction of making the circuit 
unstable. Thus a series inductance worsens stability. Of course, some 
inductances have rather large series resistances, and this may, in itself, 
tend to improve stability. 


6.2.5 Constant-Current Supply—Resistive Choke Stabilizer 


The need to get rid of the de loss in the stabilizing circuit forces one 
to consider parallel inductance-resistance networks as shown in Fig. 15. 
The constant-current supply is shown with parasitic shunt capacitance 
Cy and series resistance R,. R, must be greater than R_ — R,. for this 
scheme to work, but it is used only to stabilize the lower frequency 
range, and may therefore be put further from the diode. Between R, 
and the diode, there is additional shunt capacitance denoted by (C1. 
To stabilize the effects of Ci, a parallel inductance-resistance network 
can provide the locus Zs shown (approximately) and therefore stabilize 
the oscillator. A simple design criterion may be formulated by assuming 
R, to be very large. By choosing L/R large enough, the locus of Zp 
may be made to resemble that of the pure series resistance, shunt 
capacitance case as closely as desired. From a graphical evaluation, 
the criterion L ~ 3C,R? together with R = 2 to 5 times R_ — R,, 
gives an adequate margin of stability. For the GaAs oscillator circuit, 
an experimental model of this choke was made with R = 470 ohms, 
and L = 75uH. This would satisfy the above criterion for a capacitance 
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(a) (b) 


Fig. 15—Stability considerations for constant-current supply with resistive choke 
stability network: Section 6.2.5. 


of 100 pF’. In actual fact, it was possible to make C; = 3000 pF without 
causing instability. With C,; = 10,000 pF, the diode burned out. The 
technique used to construct the choke is illustrated in Fig. 16. The 
body of the choke is made of a ferrite material similar to that used in 
low-frequency inductors. A hole is drilled through the center and a 
470-ohm §-watt carbon resistor inserted. The winding around the 
ferrite is made of three layers of number 35 FORMVAR wire. Attached 


3 LAYERS NO.35 
SOLDER, 7 FORMVAR WIRE BRASS CAP 






~—— — Se ee Ke Ke eK YY SOLDER 


eeoererse 


Y Uj : 


YA 


\ 
470 OHM 1/8 WATT \SSOLDER 
CARBON RESISTOR 


}<-—0,100 1In— >| 


Fig. 16—Mechanical construction of the resistive choke stabilizing network. 
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to the ends are brass caps which thread into the center conductor of 
the coaxial bias line. This choke is placed at the same position indicated 
for R, in Fig. 1. A broader bandwidth might be achieved by lowering 
the resistance and inductance values, with some sacrifice of stability 
margin. 

In the particular model built, no attempt was made to insure low 
reflections at microwave frequencies from the bias circuit when the 
choke was inserted. It would be desirable to do this for a practical 
oscillator. 

With this resistive choke stabilizing network, a GaAs IMPATT 
oscillator was made to give 3.4 watts of output power with 12 percent 
efficiency at 250 mA bias current. Without any stabilization, diodes of 
this same type (breakdown capacitance ~2.5 pF’, breakdown voltage 
105 volts) consistently burned out at 80 to 100 mA. This network 
has vastly improved the stability of the bias circuit, but provides 
lossless de power transmission from the constant-current supply to the 
diode. 

Using a more sophisticated biasing network, designed upon the 
criteria developed here, but which provides a more controlled im- 
pedance locus than achievable with a simple RL choke, it has been 
possible to stabilize the bias circuit over all current ranges up through 
thermal burnout of the diode, without loss of RF or de power. 

A word of warning about the general use of chokes for stabilization : 
if the inductance is not large enough, such a choke actually worsens 
the instability and hastens burnout. 


6.2.6 Bias-Circuit Line Length 


In the previous examples it has been assumed that the terminations 
specified were all positioned exactly at the diode. In fact, this cannot 
be true and, therefore, a finite length of transmission line is necessary 
between the diode and the termination. This effect causes the termina- 
tion impedance to be transformed to a lower impedance, and increases 
the likelihood of an instability. This places a limit on the length of 
line that can be tolerated between the diode and the stabilizing net- 
work. For example, for the GaAs data of Fig. 2, and a pure stabilizing 
resistance of 500 ohms, the maximum length of line is about Imax 
=0.09\ at a frequency of about 13 MHz. But \ = c/(fve-), where 
e, is the dielectric constant of the material in the bias line. This gives 
about 2 meters for air line with e, = 1. In common microwave ab- 
sorbent materials such as Eccosorb MI 124, e, = 27 and l reduces to 
40 cm. This is for a medium Q oscillator, Q; being ~250. For a low-Q 
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oscillator (Q; ~ 25) the frequencies all will scale with 1/Q so the line 
length then reduces (in the Eccosorb material) to ~4 cm, which is 
beginning to be difficult to do. 

If the termination is frequency dependent, the effect of line length 
can be more pronounced. An example, a length / of transmission line 
terminated in a pure capacitance, is illustrated in Fig. 17. Since at 
f = fu, Xe/Z. = — 1.8 for the GaAs oscillator (see Fig. 2), the 
frequency at which X2/Z, = — 1.3 gives the maximum fy tolerable 
without initiating bias instabilities for each value of C. This value of 
fu is plotted versus C for 1 = 10, 20, and 30 centimeters, and may be 
used to derive either (7) the maximum value of shunt capacitance, 
(2t) the maximum length / for a given C and fy, or (dz) the smallest 
value of microwave loaded Q [eq. (11) ], which is tolerable before 
oscillation or instability occurs. For example, with a capacitance 


nN 
j=] 
oO 


fy IN MEGAHERTZ —> 
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Fig. 17—Depicting the effect of a length of line, 1, between the diode (having maxi- 
mum bias-circuit oscillation frequency fy) and a shunt capacity C, for the GaAs diode 
terminal impedance shown in Fig. 2. 
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C = 50pF, f. = 6 GHz, 1 = 10cm, and R_/R., = 4.3 we find fu & 47 
MHz and Q, = 232. Thus the loaded Q would have to be greater than 
232 in order that a 50-pF shunt capacitance at a distance of 10 cm 
from the diode would not cause instability in a 6-GHz GaAs oscillator 
similar to those used in these experiments. 


6.2.7 Higher-Impedance Bias Circuits 


If the bias-circuit characteristic impedance were made several times 
larger than the maximum induced value of negative resistance, then 
a matched termination at all frequencies would insure stable operation. 
Even small reflections could easily be tolerated, regardless of their 
distance from the diode. The difficulty with this is that it is very 
difficult to make good circuits with very-high-impedance lines. Tra- 
ditionally, 50 ohms is used. Any increase above that value would 
materially assist the stabilization of the bias circuit. On the other hand, 
one could also decrease the impedance of the diode by making its 
area larger. In either case, difficulty will eventually be encountered 
because of the problem of matching the diode to the microwave circuit, 
but it is probable that some gain over the present situation can be had 
by these techniques. 


Vil. SCALING LAWS 


The negative resistance of eq. (2) is obviously a function of the 
diode area and length, the maximum field strength, and the avalanche 
coefficient. It is therefore reasonable to inquire as to how R_ scales 
with different oscillator designs, material, and frequency. This can be 
worked out rather simply for the Read diode model assumed here, 
but for abrupt junction pn diodes this model fails to answer the ques- 
tion, how should the design of the diode be changed to decrease the 
value of R_, while still achieving high-power, high-efficiency micro- 
wave operation? 

For the Read diode, in the large-pulse approximation we may assume 
the negative conductance is approximately 
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Therefore, 
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6 being the normalized ac voltage amplitude and V, the breakdown 
voltage of the diode. The ratio 6 is assumed here, as elsewhere,"* to be 
invariant with frequency scaling, and with this assumption we find 
that R. for a new, optimally designed diode, as compared with R_, 
for a reference diode, is given by 


Re m Va lao m La 
a) GG)-Gye) 
Tee Mol \Vao da Mo} \ZLao 
Here Zz is defined as the ratio of the breakdown voltage to the operat- 
ing current and W-EF, has been set equal to Va. 

An interesting application of (12) is to the millimeter-wave Si 
IMPATTs designed to work at 100 GHz. For these diodes Vz = 13 
volts, and a typical current is 0.1 ampere, giving Zz = 130 ohms. For 
typical 6-GHz Si IMPATTs Vz = 105 volts, and a typical current is 
0.20 ampere, giving Zz. = 525 ohms. Therefore, 


R_(100-GHz diode) ( m ) 
ee X 0.25. 
R_,(6-GHz diode) 











Mo 


It is found that R_, is a few ohms and that m/m, < 1, so it is predicted 
that the millimeter-wave IMPATT diodes should have very small 
induced negative resistance. 


VIII. SUMMARY AND CONCLUSIONS 


In this paper we have shown that a low-frequency negative resis- 
tance is induced in the IMPATT diode by the large-signal character- 
istics of the oscillator. This low-frequency negative resistance has been 
shown to be responsible for both bias-circuit oscillations and a class 
of low-current burnouts, normally called ‘“‘tuning-induced”’ burnouts. 
The large upconversion of noise that often occurs is also caused by this 
same effect. 

The origin of the negative resistance was shown to be the combina- 
tion of the rectification property of the nonlinear ac avalanche, the 
coupling of fluctuations in the bias current to fluctuations in the 
microwave voltage amplitude by the microwave circuit oscillator 
constraint, and the fact that increases in the microwave voltage 
amplitude and de bias current generally drive the large-signal micro- 
wave admittance of the diode in opposite directions on the admittance 


300 THE BELL SYSTEM TECHNICAL JOURNAL, MARCH 1973 


plane. Calculations were made showing the correctness of this picture, 
and were applied to Si, Ge, and GaAs diodes designed for 6-GHz opera- 
tion. The agreement with experimental results is good, showing a small 
value of negative resistance for Si (typically less than 10 ohms), some- 
what larger for Ge (typically 35 to 75 ohms), and larger yet for GaAs 
(typically 85 to 150 ohms). 

A reasonably complete experimental characterization of this nega- 
tive resistance was given, showing its dependence on baseband fre- 
quency, dc bias current, and microwave loading conditions. It was 
found that, at low enough bias currents, the negative resistance was 
small for microwave loads near those which gave maximum power. 
At higher bias currents, however, the negative resistance approached 
its maximum value at or near optimum loading conditions. It was also 
found that the negative resistance had a low-pass frequency character- 
istic and an inductive reactance associated with it. The negative re- 
sistance was found to go to zero at a frequency called fm, the theo- 
retical maximum frequency of oscillation, and fw was found experi- 
mentally to be roughly inversely proportional to the microwave cir- 
cuit’s loaded Q factor. It was found that, for large-signal bias-circuit 
oscillations, the negative resistance decreased, and the reactance 
remained approximately constant. This fact indicates open-circuit 
stability, and is very important in the stabilization of the bias circuit. 

An analysis was performed and an equivalent circuit was derived, 
showing that the inductive reactance observed experimentally is best 
represented by a negative capacitance, —1/yR_, in parallel with the 
negative resistance, — R_, all in series with the space-charge resistance, 
R,-. This circuit predicts several aspects of the small-signal behavior of 
the induced negative resistance, and allows the stabilization tech- 
niques to be designed in a straightforward manner. The agreement of 
this equivalent circuit with the experimental results is good and it 
shows in detail how the interaction of the microwave circuit and diode 
admittances affect the low-frequency negative resistance. In particular, 
it is shown that the maximum frequency of oscillation, fy, is given by 


R_ 
fu = = Pe 
Qr Rec 
where f, is the microwave frequency of oscillation, Qz is the loaded 
Q of the microwave cavity, R,. is the space-charge resistance of the 
diode, and R_ — R,, is the low-frequency asymptote of the negative 


resistance. Smith chart plots of the negative of the diode’s small- 
signal terminal impedance, (—)Z;, are given for various values of 
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R_/R;_ and normalized frequency. The shape of these loci is found to 
vary somewhat, though not drastically, with the choice of impedance 
normalization used. When normalized to 50 ohms (which is greater 
than R,.), the contours always appear somewhat as shown in [ig. 2. 

Stabilization techniques were discussed, beginning with the develop- 
ment of a simple means of judging the stability of a given oscillator 
and bias circuit. Several examples of stable and unstable conditions 
were given including: 


(z) a calculation of the maximum equivalent lumped shunt capaci- 
tance that can still remain stable, 

(22) discussion of the effects of line length in the bias circuit which 
indicated the desirability of placing any stabilization network 
as close to the diode as possible, but at the same time indicating 
that it need not be exactly at the diode, 

(i2t) a stabilizing network consisting of a pure series resistance, R,, 
in the bias circuit with R, > R.— R,. (it is required that 
there be no significant reactance or line length between the 
diode and R;), and 

(tv) an alternative stabilizing network consisting of a parallel in- 
ductance-resistance choke combination that does not dissipate 
de power. 


Using the principles of stabilization outlined here, it has been possible 
to stabilize the bias circuit of IMPATT oscillators over the full range 
of bias currents up to the diode’s thermal burnout limit. 

As a final topic, the frequency scaling of R_ was considered, and it 
was shown as an example that the 100-GHz Si IMPATTs should have 
very small induced negative resistances, and therefore the instabilities 
encountered at lower frequencies and in Ge and GaAs should not be 
a serious problem for the millimeter-wave Si diodes. 
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APPENDIX 


It is the purpose of this Appendix to derive the equivalent circuit 
of the low-frequency diode impedance of an IMPATT diode oscillator. 
To do this, the Read diode model is assumed to describe the diode and 
a quasi-static technique is assumed to describe its interaction with the 
microwave circuit. Rewriting eqs. (8) and (9) of Kurokawa’s work" 
to apply to admittance, and dropping the noise voltage terms, we may 


write Av 
[G(w) — 9]B'(e) — (Be) + 616") + 1¥'@) |? =F re (18) 








d 
[G(w) — g]G"(w) + LB() + 6B’) + [¥"(w) |? = 0. (14) 


Here, Y(w) = G(w) + jB(w) is the microwave circuit admittance, 
y = —g + jb is the device admittance, V. and ¢ are the microwave 
voltage amplitude and phase, and are assumed to be slowly varying 
functions of time, and the prime on G, B, and Y denotes differentiation 
with respect to frequency. Equations (13) and (14) describe the change 
in time of the amplitude and phase of a negative conductance oscillator. 
In steady-state oscillation, dV./dt = 0 and dy/dt = 0, implying 


G(w) ps g(V.,f) = 0, (15) 
and 
Bw) + b(V,,J) = 0, (16) 


which determine the steady-state oscillation amplitude A, and fre- 
quency w.. We have assumed that g and b are independent of frequency, 
which only implies that the analysis is restricted to situations in which 
Y(w) varies much more rapidly than y. 

Let us now presume that the steady-state amplitude and bias cur- 
rent are perturbed by a small amount, that is 


Ig=I,+Al 
and 
V.=A,+A6 


where AI «J, and 6 <1. Then, linearizing the variation of g and b 
with both V, and Ja gives 


g=G-—sG6+ k,Al 
and 
b= —B+rGi — k,AI 
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where 
A. og A. ob 
een Ge Oi GeV, 
(17) 
2 og eee ab 
A ape ° al 


are the RF voltage and de current linearization parameters of the 
diode’s admittance, evaluated at the steady-state large-signal oscil- 
lation point. We may then write 


dé 
[s@6 — k,AI]B’ — [rG6 — bya] + |¥"/?— = 0. 


Collecting terms, 


dé 
ate vO Cay. (18) 
at 
where 
sGB’ — rGG’ 19) 
"G2 BP 
and 
B'k, — G’ky 
= (20) 
G2 ++ BY? 
Introducing 
G'/B’ = tan 6, 
r/s = tan 6a, 
and 
6= 6. + 6a+ 7/2, 
we obtain 
Vr2 + 3? 
y = ——— Gsin 8, (21) 
oY 
Ow 








which is eq. (7) of the text. 
Equation (19) for Af becomes 





: Pane 
I eee 
ay 


Ow 


sin (x)Al. (22) 








The meaning of the angles @., 62, and 6 can be seen in Fig. 18 where 
it is seen that @ is the total angle of intersection between the negative 
of the circuit admittance curve, and the tangent to the large-signal 
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0 





Fig. 18—Defining the angles 6,, 64, and @ in the microwave large-signal admittance 
plane of the diode. 


device variation with voltage amplitude, at the steady-state operating 
point. The angle x is defined by 


Tv 
X= Or — 8 
with 6, defined as above and 6; defined by 


k 
tan Oy = ee 
9g 


The graphical interpretation of x is given by Fig. 19, and it is seen 
to be the angle between the negative of the circuit admittance and the 
direction of the variation of the large-signal device admittance with de 
current, also at the steady-state operating point. 
Assuming an exp(jwt) time dependence for 6 and AI, that is 6 
= Real (de%*), and AJ = Real (i,,e%*), eq. (17) gives 
Pa ee 
y(1 + jw/y) 
or 
Veo + ky sin x 1 im 


Vr? 4 §2 Sin 6G (1 + jo/y) 
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Fig. 19—Defining the angles @7, @., and x in the microwave large-signal admittance 
plane of the diode. 
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~ Agsin ol dla {(1 + jw/y) 








We may therefore interpret y as the modulation frequency at which 
the microwave power in the sideband (~|6|*) is one-half of its low- 
frequency response, that is, the common 3-dB modulation half-band- 
width of the oscillator. 

The rectification equation, eq. (1) of the text, is now assumed to 
be true as V, varies slowly with time. This gives, for the voltage at 
frequency w (to first order in 4), 





~ m 
Va — Reclm nae 2V.6 
4W 


or 

dvaj 1 
fer (23) 
AEs Guy) 


Vin R m sinx 
= SC 6 0 Se 
im 4We,sin 6 
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We may then write 


i 


Z: = Ree — —————_ 
1+ jw/y 


(24) 


where R_ and y are defined by eqs. (6) and (7) of the text. Z, is the 
terminal impedance of the diode at baseband frequencies, and if 
R_ > R,-, then a net negative resistance exists at low frequencies. The 
equivalent circuit for Z,; is shown in Fig. 7 and its properties are dis- 
cussed in Section V of the text. 


REFERENCES 

1. Clorfeine, A. 8., and Hughes, R. D., “Induced de Negative Resistance in Ava- 
lanche Diodes,” Proc. IEEE (Letters), 57, No. 5 (May 1969), pp. 841-842. 

2. Olson, H. M., Jr., “Negative Resistance Diode Coaxial Oscillator with Resistive 
Spurious Frequency Suppressor,’ U. 8. Patent No. 3,621,463, November 16, 
1971. 

3. ‘“Microwave Power Generation and Amplification using IMPATT Diodes,”’ 
Hewlett-Packard Application Note 935, Hewlett-Packard Corp., Palo Alto, 
California. 

4. Read, W. T., Jr., ““A Proposed High-Frequency, Negative Resistance Diode,’’ 


11. 
12. 


13. 
14. 


B.S.T.J., 37, No. 3 (March 1958), pp. 401-446. 


. Evans, W. J., Scharfetter, D. L., Johnston, R. L., and Key, P. L., ‘Tuning 


Initiated Failure in Avalanche Diodes,” J. Appl. Phys., 42, No. 2 (February 
1971), pp. 799-803. 

Blue, J. L., “Approximate Large-Signal Analysis of IMPATT Oscillators,” 
B.8.T.J., 48, No. 2 (February 1969), pp. 383-396. 


. Sze, S. M., Physics of Semiconductor Devices, New York: John Wiley and Sons, 


1969. 


. Gewartowski, J. W., and Morris, J. E., “Active IMPATT Diode Parameters 


Obtained by Computer Reduction of Experimental Data,’ IEEE Trans. 
Microwave Theory and Techniques, MTT-18, No. 3 (March 1970), pp. 
157-161. 


. Magalhaes, F. M., and Kurokawa, K., ‘‘A Single-Tuned Oscillator for IMPATT 
10. 


Characterizations,’’ Proc. IEEE (Letters), 48, No. 5 (May 1970), pp. 831-832. 

Gupta, M. S., ‘“Noise in Avalanche Transit-Time Devices,’’ Proc. IEEE, 59, 
No. 12 (December 1971), pp. 1674-1687. 

Kurokawa, K., ‘‘SSome Basic Characteristics of Broadband Negative Resistance 
Oscillator Circuits,’ B.S.T.J., 48, No. 6 (July-August 1969), pp. 1937-1955. 

Scharfetter, D. L., and Gummel, H. K., “Large-Signal Analysis of a Silicon Read 
Diode Oscillator,’ IEEE Trans. Electron Devices, HD-16, No. 1 (January 
1969), pp. 64-77. 

Kurokawa, K., “Noise in Synchronized Oscillators,’ IEEE Trans. Microwave 
Theory and Techniques, MTT-16, No. 4 (April 1968), pp. 234-240. 

Scharfetter, D. L., ‘‘Power-Impedance-Frequency Limitations of IMPATT 
Oscillators Calculated from a Scaling Approximation,’’ IEEE Trans. Electron 
Devices, HD-18, No. 8 (August 1971), pp. 536-548. 


Copyright © 1973 American Telephone and Telegraph Company 
THE BeuL SysTEM TECHNICAL JOURNAL 
Vol. 52, No. 3, March, 1973 
Printed in U.S.A, 


Dynamic Data Reallocation in 


Bubble Memories 


By P. I. BONYHARD and T. J. NELSON 
(Manuscript received October 10, 1972) 


Bubble technology offers several operations that have no equivalents in 
technologies based on magnetic recording. Examples of such operations 
are: transfer, reversal of the direction of propagation, and opening and 
closing of gaps in the data stream. This paper* shows how such operations 
can be used to dynamically reallocate data in the bubble memory, causing 
it to become an integrated memory hierarchy. A considerable tmprovement 
in performance results. A model is presented which relates the bubble 
memory with dynamic reallocation to stack processing, a technique used 
in the evaluation of memory hierarchies. With the aid of this model it 
becomes possible to calculate the performance of the bubble memory using 
published data derived from the traces of selected typical programs. 
Memory design ts optimized for the execution of such programs. Design 
parameters are proposed for a 2-Mb bubble memory with 128 detectors 
which, in the execution of the type of program for which data were available, 
requires an average of only 8.8 shifts for access and an average of 12.1 
shifts per memory cycle. If bubbles are propagated at a rate of 1 MHz, the 
average access and cycle times for this memory become 8.8 us and 12.1 us, 
respectively. Such performance, in conjunction with the low cost per bit 
offered by bubble technology, 1s expected to have a major impact. The per- 
formance of this memory, when operated in conjunction with a faster 
buffer, is also calculated. The use of a 64-kb buffer 1s shown to reduce the 
average number of shifts for access to 1.05, and the average number of 
shifts per cycle to 1.9. 


I. INTRODUCTION 


One major area of application of magnetic bubble technology is 
mass memory.! The potential advantages of bubbles over disk files 
and drums are: shorter access time, lower cost and power dissipation, 

* The contents of this paper were presented at the 1972 Intermag Conference by 
P. I. Bonyhard. 
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reduced volume, and a higher degree of modularity.? Moreover, there 
are several operations which can be performed on the information in a 
bubble memory which cannot be performed by conventional magnetic 
recording. The purpose of this paper is to show that a magnetic bubble 
memory can function as a hierarchy, given certain operations which 
are not difficult to realize. 

The operations relevant to this paper are: transfer, instantaneous 
reversal of the direction of propagation, and removal of bubbles from 
a propagate channel while closing the gap left by them or creation of 
a gap while inserting new bubble information into the channel. Transfer 
of bubbles from one propagate channel to another already plays an 
important role in the design of bubble mass memories.! By a process 
of (2) propagating forward n cycles, (227) removing bits from the channel 
and closing the gap, (77) propagating backward n cycles, and (iv) 
opening a gap and reinserting the bubbles, a permutation of the 
information in the channel is effected. It will be shown how a simple 
algorithm using one permutation element per propagate channel causes 
the memory to function as an integrated hierarchy. We also show how 
this dynamic reallocation of data can be combined with the major- 
minor loop organization! so as to optimize memory cost performance. 
It is concluded that the improvement is so substantial as to have a 
major impact in the computer industry. 


II. DYNAMIC DATA REALLOCATION IN CLOSED LOOP SHIFT REGISTERS 


Consider an assembly of simple, closed loop shift registers as shown 
in Fig. la. Information propagates in all registers synchronously under 
the influence of a common rotating drive field. One page of information 
is stored along a horizontal line, a particular page being formed by 
the black dots in Fig. la. It is inherent in the dynamic scheme that 
each page must carry its own reference address and some of the bits 
of the page are devoted to this purpose. Thus, extra loops, totaling 
loge (number of pages), have to be provided. 

Upon request for a specified page, information is shifted clockwise 
until the appropriate group of detectors, the position of which is 
marked D in Fig. 1b, detects the right address. This takes, say, x cycles 
of propagation. Now the page can be read or rewritten. The memory 
next is reset by shifting x cycles counterclockwise. Using the circuit 
arrangement of Fig. 1b, the addressed page will remain arrested at 
the detectors while all other pages move back. A T-bar realization 
of Fig. 1b is shown in Fig. 2. This design has been operated quasi- 
statically to demonstrate feasibility. 
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(a) (b) 
Fig. 1—-Assembly of shift register loops for dynamic address reallocation. 
Now let the physical page locations be numbered 1, 2, 3, --: , n, 


according to their distance from the detectors. A given page will reside 
in a location, the number of which is equal to the number of requests 
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Fig. 2—Magnetic circuit used to realize the function of Fig. 1b. 
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that have been made, since the page was last requested, for pages 
originally in locations with higher numbers. Pages never requested will 
lie in the locations with the highest numbers. Thus recently used 
pages are near the ‘‘top,’’ whereas less recently used pages are 
farther down. This replacement algorithm has been discussed before 
in the literature and has been named “‘the stack.’’ 

It will now be shown how the average number of shifts necessary 
to reach an addressed page in the bubble stack can be calculated. 
Consider the stack to be arbitrarily divided into two parts, one part 
consisting of pages in locations 1, 2, --- , ki, and the other part of 
pages in locations k; + 1, k; + 2, --- , n. It should be recognized that 
the first k; page locations can be thought of as forming a “‘buffer’’ and 


TaBLeE I—Hir Rario Data (REPRODUCED FROM REF. 4) 





Buff 
Size ; Page Size (Bits) 


(Bits) Classes} 128 256 512 1k 2k 4k 8k 16k* 32k* 
8k 1 0.894 0.915 0.928 0.924 0.884 0.792 0.495 


4 0.895 0.916 0.921 0.904 0.791 — == 
16 0.891 0.903 0.860 — — = — 


Pe 
Lied 


64 0.857 — ae — = _ 
16k 1 0.931 0.949 0.957 0.958 0.950 0.912 0.824 0.56 — 
4 0.931 0.948 0.954 0.955 0.933 0.808 — — tre 
16 0.930 0.943 0.943 0.909 — a — a tee 
64 0.921 0.913 — _ = — = ee 
32k 1 0.951 0.969 0.973 0.978 0.977 0.966 0.939 0.86 0.64 
4 0.955 0.969 0.973 0.977 0.974 0.951 0.884 — — 
16 0.955 0.968 0.972 0.970 0.933 — — =. seme 
64 0.955 0.963 0.948 — as _ _ eee “Ga 
256 0.934 — — _ _ — — —- — 
64k 1 0.977 0.986 0.988 0.985 0.987 0.987 0.984 
4 0.981 0.986 0.988 0.986 0.987 0.985 0.965 


16 0.981 0.985 0.988 0.987 0.983 0.954 
64 0.979 0.984 0.985 0.974 — — 
256 0.974 0.971 — = = = 


128k 1 0.985 0.993 0.994 0.996 0.993 0.992 0.994 
4 0.990 0.993 0.994 0.996 0.994 0.992 0.993 

16 0.990 0.994 0.995 0.997 0.995 0.991 0.957 

64 0.990 0.994 0.995 0.995 0.985 — = 

256 0.989 0.992 0.986 — = = rae 


256k 1 0.989 0.996 0.997 0.997 0.999 0.994 0.997 
16 0.994 0.996 0.997 0.998 0.998 0.996 0.997 
32 0.994 0.996 0.998 0.998 0.998 0.997 0.997 
64 0.994 0.996 0.998 0.998 0.998 0.988 — 
256 0.994 0.996 0.997 0.992 — = —_ 


DS a 
tA ar Ws 


i a 
Re eee Ee 
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the last n — k; pages a ‘‘memory”’ in the terminology conventionally 
used for two-level memory hierarchies.?# Whenever a page is brought 
from the memory to the buffer, the page currently in the kth position 
moves from the buffer to the memory. Clearly, this page is the least 
recently used page in the buffer, so that the replacement algorithm in 
this two-level hierarchy is “least recently used’ (LRU). Hit ratios, 
that is, fractions of all requests that can be satisfied from the buffer 
without reference to the memory, can be found in literature. A par- 
ticularly useful set of hit ratio data is given in Table II of Ref. 4 and 
is reproduced as Table I of this paper. The columns marked with 
asterisks contain entries obtained by graphical extrapolation from the 
original data. 

In a two-level hierarchy, derived by cutting the bubble stack at 
the k;th location, the number of bits per page is equal to the number 
of bubble loops, say, ¢. The number of bits in the buffer is then ¢k,. If 
the corresponding hit ratio is h;, then the average number of shifts 
necessary per request might be estimated as 


- k-1 k; -—1 
se yp = ay, 





This is the average number of shifts in the buffer times the probability 
of a hit, plus the average number of shifts to the memory times the 
probability of a miss. 
The bubble stack, however, can be cut at any location. A better 
estimate is obtained by dividing the stack into m > 2 levels 
"ke +kei-—1 


S = 2 . (hi — hi-1) (1) 





where k, = 0, km = n, and hz = 1. The best approximation to 8S is 
achieved when every value of k between 1 and n is taken into account: 


SD (k — 1)(hi — hin). 
k=1 
This result is, of course, self-evident, as k — 1 is the distance of the 
kth location and h;, — hx_-1 1s the probability of hitting the kth location. 
S has been calculated on the basis of eq. (1) using the entries that 
can be found in Table I. It is an overestimate because it divides the 
probability equally among the levels between entries. Actually, the 
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Fig. 3—Schematic representation of the operation of a major-minor loop organized 
memory with dynamic data reallocation in the major loops. 


probability must be monotonically decreasing. The results are tabu- 
lated below: 


Number of loops 128 256 512 1K 2K 4K 8K 
Minimum number of bits per loop 2K 1K = 512 256 128 64 32 
Average number of shifts 57.8 25.9 12.13 567 3.08 1.44 1.05 








The minimum number of bits per loop listed in the above tabulation is 
based on the information kindly provided by the author of Ref. 4 to 
the effect that the average program size used in deriving the hit 
ratio data was 256 kb. Thus closed loop shift registers can access data 
at least ten times faster with dynamic reallocation. A magnetic bubble 
memory with this feature is, in fact, an integrated hierarchy. 


III. DYNAMIC DATA REALLOCATION IN MAJOR-MINOR LOOP ORGANIZED 
MEMORY PLANES 


Dynamic data reallocation in closed loop shift registers does not 
appear to be a very attractive bubble mass memory design. Each loop 
must have its own detector, so that either the number of detectors or 
the average number of shifts is high. The average number of shifts for 
a given number of bits per detector does not compare too favorably 
with the major-minor loop organization.! Dynamic reallocation may, 
however, be combined with major-minor loop organization to form an 
extremely attractive design. Dynamic reallocation can be performed in 
the major loops only, in the minor loops only, or in both. The last of 
these three options is probably too complex to be of practical use. A 
design based on the first option is presented below. The second option 
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may be a valid alternative but, at least in the authors’ opinion, it is 
less attractive. 

Figure 3 illustrates schematically how the system operates. The 
information that may occupy the assembly of all the major loops 
concurrently forms a class. One of the classes is selected by shifting 
information in the minor loops until the required class is aligned for 
transfer, and then transferring it into the major loops. Two bits from 
each minor loop are transferred and the major loop is filled completely. 
Consequently, two bits of each class are stored in each minor loop in 
each memory plane. A plane consists of one major loop and all the 
minor loops associated with it, as well as one detector and one con- 
trolled bubble generator operating on the major loop. Once the right 
class resides in the major loops, the major loops operate exactly as the 
assembly of loops described in the previous section. The net result 
is that page locations are dynamically reallocated within each class, 
but pages are never permitted to cross class boundaries. Hit ratios 
that correspond to such a multiclass system are given in Ref. 4, and 
also appear in Table I. 

It should be recognized that a page in this system consists of all 
those bits that may concurrently be detected. Thus there is one bit of 
each page in each plane. Also, the transfer mechanism considered here 
is of the “‘conductor”’ variety,! so that reversing the sense of propaga- 
tion may be freely used to accomplish dynamic address reallocation 
in the major loops. As this consists of an equal number of forward and 
reverse shifts, the gaps left in the minor loops are correctly aligned at 
the time the class is to be transferred back. 

Reaching a page in the memory is accomplished in two steps. First 
the right class is positioned and transferred into the major loops, and 
then the page is brought to the read/write port. The latter step has 
already been discussed in the previous section, but the former one, 
class swapping, needs some elaboration. Following the approach taken 
in Ref. 4, it is assumed that several programs are concurrently resident 
in the memory. More specifically, the memory size is chosen to be 
2 Mb, whereas the average program size is 256 kb, so that an average 
of 8 programs may share the memory. It is assumed that each program 
resides in some number of contiguous classes which will be called 
active classes for the program currently being executed. The class 
occupying the major loops at any given time can be looked upon as a 
single-page buffer of a single-class two-level hierarchy, with all the 
other classes forming the rest of the pages in the memory. Of course, 
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for a single-page buffer the replacement algorithm is trivial. However, 
as can be seen in Table J, the hit ratios associated with such single- 
page buffers are still fairly high. Consequently, there is a good chance 
that the next request is addressed to the class already in the major 
loops, provided that the addressed class is not transferred back to the 
minor loops until it is known which class is addressed next. Therefore, 
it makes good sense to operate in this manner. 

The average number of shifts necessary to bring out the new class, 
provided that a ‘‘class-page’’ failure has occurred, will now be calcu- 
lated. It is assumed that the particular program executed at the time 
resides in m contiguous classes and the 7th class was addressed last. 
In the absence of any further information, it appears to be reasonable 
to assume that each other class is equally likely to be addressed next. 
Consequently, the average number of shifts is 

22 “(55+ : i) (4 mar a+1)(m t) 


j=1 j=l m 


It has to be remembered that there are two bits per class in each minor 
loop. The average S must be further averaged as 7 assumes all values 
between 1 and m, so that 


2m? — 1 


ae 





sa — als (ij ~ 1)t = 
i=1 m? t=1 

To this two further shifts must be added to accomplish in and out 
transfer of two bits per minor loop. The total average number of class 
swapping shifts per page request is 





2m? — 1 
Su = (1 hed (5 . +2) (2) 


where h,, is the class hit ratio discussed earlier in this section, and is 
the number of active classes. 

To the class swapping shifts, given by eq. (2), one must add the 
average number of shifts necessary to reach a page within the selected 
class to get the total average number of access shifts per request to 
the memory. The average number of shifts within the selected class is 
given by eq. (1) with the hit ratios taken for the appropriate number 
of classes. The latter quantity must be doubled when calculating the 
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average number of shifts per memory cycle. The results are given 
below for three alternative designs labeled A, B, and C. 








A B C 
Bits per memory 2M 2M 2M 
Planes per memory 128 128 128 
Bits per plane 16k 16k 16k 
Classes 64 128 256 
Bits per minor loop 128 256 512 
Minor loops per plane 128 64 32 
Bits per major loop 256 128 64 
Active classes 8 16 32 
Active bits per minor loop 16 32 64 
Class hit ratio = 0.64 = 0.56 0.495 
Class swapping shifts 7.25 12.63 23.3 
Class swapping shifts per request 2.61 5.56 11.8 
Page search shifts 7.00 3.26 1.40 
Total shifts per access 9.61 8.82 13.2 
Total shifts per cycle 16.61 12.08 14.6 








Details of the page search shift calculations are given in Table II. 
The memory size was chosen to be consistent with Ref. 4, where the 


Taste II—Ca.LcuLaTION OF THE AVERAGE NUMBER OF PAGE 
SEARCH SHIFTS FOR THREE DESIGNS 





kis + hin 1 3 (hi — hit) X 
Design ki 2 h; hi — hie (ki + kia — 1) 

8 3.5 0.894 0.894 3.13 
16 11.5 0.931 0.037 0.43 
A 32 23.5 0.955 0.024 0.56 
64 47.5 0.981 0.026 1.24 
128 95.5 0.990 0.009 0.86 
256 191.5 0.994 0.004 0.77 
6.99 
4 1.5 0.891 0.891 1.34 
8 5.5 0.930 0.039 0.21 
16 11.5 0.955 0.025 0.29 
B 32 23.5 0.981 0.026 0.61 
64 47.5 0.990 0.009 0.43 
128 95.5 0.994 0.004 0.38 
3.26 

2 0.5 0.880 0.880 0.440 

4 2.5 0.934 0.054 0.1385 

C 8 5.5 0.955 0.021 0.115 

16 11.5 0.980 0.025 0.287 

32 23.5 0.990 0.010 0.235 

64 47.5 0.994 0.004 0.190 


1.402 
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data came from, and the plane size was chosen as the one that would 
provide an economical design on the basis of current cost estimates. 

Each of the three designs has an overhead associated with it, because 
of the extra planes needed to store the address tags of the reallocated 
pages. This overhead for designs A, B, and C is 6.3 percent, 5.5 percent, 
and 4.7 percent, respectively. 


IV. THE USE OF A BUFFER’ TO IMPROVE PERFORMANCE 


The best design outlined in the preceding section gave an average 
cycle time of about 12.1 shifts or, assuming 1-MHz operation, 12.1 us. 
The question arises of how this figure may be further improved by the 
use of a smaller submicrosecond cycle time core or integrated circuit 
buffer. The answer can be readily calculated, provided that the buffer 
is divided into as many classes as there are active classes in the memory, 
and that the page size in the buffer is the same as in the memory. Now, 
if there are kg page frames per class in the buffer, then the pages 
currently residing in location 1, 2, --- , kg of all active classes will 
also appear in the buffer. This statement is not quite correct if the 
buffer stores only read (fetch) data, whereas write (store) data are 
sent from the central processor directly to the memory, but the differ- 
ence is probably negligible. 

It is, of course, also necessary to uniquely assign all nonactive 
memory classes to buffer classes, so that program swapping can_take 
place. This can be done as follows. In terms of Design B, let the 7 most 
significant bits of the 14-bit page address be the class address. As 
each program occupies contiguous classes, the 4 least significant bits 
of the class address refer to classes in the same program for a typical 
program size. Thus the Ist, 2nd, 38rd, and 4th least-significant bits of 
the memory page address should be considered to be the buffer class 
address. 

The average number of shifts necessary in the memory per request 
can now be calculated as follows. If the buffer hit ratio is hz, then hes 
in eq. (2) should be replaced by hz, whereas all contributions to the 
page search shifts, calculated in Table II, by values of k, S$ kz should 
be neglected. The results are tabulated below, still for Design B. 


Buffer size (bits) 8k 16k 32k 64k 
Buffer bit ratio 0.891 0.930 0.955 - 0.981 
Pages per class in buffer 4 8 16 32 
Class swapping shifts 1.37 0.88 0.56 0.24 
Page search shifts 1.92 1.71 1.42 0.81 
Total access shifts 3.29 2.59 1.98 1.05 
Total shifts per cycle 5.21 4.30 3.40 1.86 
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V. DISCUSSION AND CONCLUSION 


The availability of memories costing approximately the same or 
less per bit as a disk file or drum with cycle times of about 10 us is 
likely to have a major impact on the computer industry. This appears 
to be the answer to the system designer’s old dream of a memory at 
core speeds and disk costs. It looks like we can fulfill this dream with 
bubble memories, provided that we can shift bubbles at 1-MHz rates, 
which appears to be a reasonable objective. It may be added that 
the organizations discussed here should be also applicable to other 
forms of serial memory technologies, MOS registers, charge coupled 
registers, ete. 

The major applications for bubble memories with dynamic address 
reallocation are likely to be in mini- and midi-computers and in other 
systems with relatively less powerful central processors. The larger 
systems and those with faster processors will probably find the buffered 
configuration more attractive. For instance, one may use a 64-kb, 
0.5-us access time, l-us cycle time buffer to give an average access 
time of about 1.5 us, and an average cycle time of about 3.0 us for 
the 2-Mb bubble memory. In many systems this would be vastly 
preferable to an all-integrated-circuit memory costing much more for 
the same capacity, even if the latter memory has a cycle time of 0.1 us. 

Further work should be directed towards the assessment of the 
applicability of these techniques to electronic telephone switching 
systems. 
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Restoring the Orthogonality of Two 
Polarizations in Radio 
Communication 


Systems, II 


By T. S. CHU 
(Manuscript received October 3, 1972) 


The Poincaré sphere has been applied to the analysis of orthogonalizing 
two polarization ellipses by a differential phase shifter and a differential 
attenuator. The condition of minimum differential attenuation for removing 
a given amount of nonorthogonality is determined. A previously reported 
transformation via two nonorthogonal linear polarizations should be used 
for two slender ellipses.1 Another transformation via two oppositely 
rotating ellipses having parallel axes and equal axial ratios should be 
used for two fat ellipses. System applications of the transformations are 
discussed. 


I. INTRODUCTION 


A method of recovering the orthogonality of two polarizations in a 
radio communication system was presented recently.1 Two arbitrary 
polarization ellipses are first transformed simultaneously into two 
nonorthogonal linear polarizations by a differential phase shifter, and 
then the nonorthogonality is removed by a differential attenuator. 
It is of interest to ask whether that transformation is optimum for 
system applications. 

Before proceeding further, we will define the optimum transforma- 
tion. Clearly, it is desirable to minimize the differential attenuation 
that we must introduce to correct for nonorthogonality. However, in 
order to achieve maximum bandwidth, we also should minimize the 
differential phase shift even if we assume, as we do here, that the 
polarization characteristics of components in the system are not very 
frequency sensitive over the operating bandwidth. This assumption is, 
for example, expected to be valid for any depolarization in the main 
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beam of reflector-type antennas provided there is no polarization dis- 
tortion in the feed radiation. If necessary, we can apply polarization 
correction separately to each of the subbands within a wide band. 

In a practical dual-polarization radio system, the two orthogonal 
polarizations feeding the transmitting antenna are either linear or 
circular. Signal generators and detectors are always linearly polarized ; 
however, conversion to circular polarization for radio transmission is 
sometimes made to avoid effects such as Faraday rotation. If the 
transmission medium and the radiating systems introduce only mod- 
erate polarization distortion, the dual polarization signals appear as 
two slender or two fat polarization ellipses at the receiving terminal, 
depending on whether linear or circular polarizations are being used. 
This classification into two types of elliptical shapes suggests an opti- 
mum transformation for each type. However, the validity of the trans- 
formations presented in the next section is independent of the shapes 
of the ellipses. A practical design of adjustable differential phase- 
shifters and attenuators has been suggested by E. A. Ohm.? 

The following analysis will be presented with the aid of the Poincaré 
sphere. This geometrical representation of the polarization of a plane 
electromagnetic wave was introduced to radio engineers by G. A. 
Deschamps.’ For convenience of the reader, a summary of the Poincaré 
spherical representation is given in the Appendix. 


II. ANALYSIS 
2.1 Minimum Differential Attenuation 


Let us first find the condition for minimum differential attenuation 
required for removing a given amount of nonorthogonality between 
two polarizations. Each polarization is represented by a point on the 
Poincaré sphere. If the great circle arc connecting two points on the 
sphere is a semicircle, then the two polarizations are orthogonal to 
each other. The degree of nonorthogonality between two polarizations 
can be measured by the deviation of the great circle arc from a semi- 
circle. Let two nonorthogonal elliptically polarized waves be repre- 
sented by points M, and M, on Poincaré sphere as shown in Fig. la. 
The great circle arc connecting M, and M, intersects the equator at C. 
The longitudes of Mi and M, with respect to C are twice the orientation 
angles of the two polarization ellipses with respect to the X-axis of a 
set of X-Y coordinates. This relationship is sketched in Fig. 1b. Since 

oN 


M;, and M:, are nonorthogonal to each other, M,C + M.C <7. Orthog- 
oy 


onalization is accomplished by stretching the great circle arc MiM, 
to the value of 7. 
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SS 
— 


Fig. 1—Orthogonalization by differential attenuation. 


Lo, oN 
tan3M,C tan3M.2C wes 
<_< = ——;_ is imposed 
tan$MiC tan3M3C 

CAN c™ 
along the X-axis such that MiC + M3C = z, then the two polariza- 
oN oN 
tions M, and M, will be orthogonalized. Since tan $MiC tan 4M3C = 1, 
the required differential attenuation is minimized by maximizing 


If a differential attenuation of 


a aT aoe 
‘Vian 3M,Ctan4M.C with a proper differential phase shift, which 
o™ o™ cr 
implies the constraint that MiC + M.C = M,M:2. In this way, the 
oN 
minimum differential attenuation needed is found to be tan 7M,M, 


a oN o™ 
when M,C = M.C. For the given degree of nonorthogonality, 7-M.Ms,, 
the orthogonalization can be performed by the minimum differential 

a 


attenuation tan {MiMz, only if the two elliptical polarizations have 
the same axial ratio. Thus the two arbitrary elliptical polarizations 
must first be transformed by a differential phase shifter to allow a 
minimum differential attenuation. 

Furthermore, one wishes to obtain two orthogonal linear or circular 
polarizations immediately following orthogonalization. This require- 
ment determines the prerequisite condition that two nonorthogonal 
elliptically polarized waves should first be transformed into two non- 
orthogonal linear polarizations or two oppositely rotating elliptical 
polarizations having parallel major and minor axes and equal axial 
ratios, both of which are limiting cases of two elliptical polarizations 
with the same axial ratio. 
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(a) (b) M;B=BM’, 


Fig. 2—Simultaneous transformation of two elliptical polarizations into two linear 
polarizations. 


2.2 Two Slender Ellipses 


If the transmitting antenna is fed by two orthogonal linear polari- 
zations, moderate polarization distortion in a radio communication 
system will produce two slender polarization ellipses at the receiving 
terminal. At that point, the orthogonalization should begin by a 
simultaneous transformation of the two polarization ellipses into two 
linear polarizations. This transformation has been obtained previously ;! 
however, it will be described here in terms of the Poincaré sphere. 

Let two nonorthogonal elliptic polarizations be represented by two 
points, M; and M;z in Fig. 2a. The intersection C of the great circle 
arc with the equator designates a set of coordinate axes X-Y which is 
shown in Fig. 1b. The polarization ratios of the two polarization ellipses 
in terms of these X-Y coordinates will be related by tan ¢1 = tan ¢2. 
Replacing each side of this equation by eq. (11) yields: 


tan 2a; esc 28; = tan 2a» ese 2Bo. (1) 


Substituting B. = 6, — 6 into eq. (1), one obtains the solution for the 
orientation of the X-Y axes 





tan 2a 
cos 26 — 
B to tan 2a1 0 < B < is (2) 
= -—cot! |———_——_——___ , —. 
‘9 sin 20 9 


oN 
The are MM, can be rotated onto the equator by applying the follow- 
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ing differential phase delay to the components in the Y direction 
@ = tan7! [tan 2a; ese 26; ]. (3) 
Now the angle » between the two linear polarizations represented 
oN oN 


by Mj and M3 in Fig. 2b is half of the sum of the ares CM; and CM, 
vent (4) 


(1 + tan? a;) — (1 — tan? a;) cos 26; 
(1 + tan? a;) + (1 — tan? a;) cos 28; 
using eq. (10). This angle y may be changed to 7/2 if a differential 
attenuation of tan (W/2) is imposed on the components in the direction 
B bisecting the two linear polarizations. This direction will be oriented 
at an angle 


where yy; = tan7! is obtained 





x = 3(v1 — ¥2) (5) 


with respect to the X direction of the coordinates in Fig. Lb. 

The above equations can be easily identified with those in Ref. 1. 
One notes that the transformations are valid for the two polarizations 
located on the same side as well as on opposite sides of the equator 
of the sphere. 


2.3 Two Fat Ellipses 


If the transmitting antenna is fed by two orthogonal circular polari- 
zations, two fat polarization ellipses will appear at the receiving ter- 
minal if the radio communication system is moderately contaminated 
by polarization distortion. Here the orthogonalization should be 
started by the simultaneous transformation of two given ellipses into 
two oppositely rotating ellipses having parallel major and minor axes 
and equal axial ratios. One looks for an X’-Y’ cartesian coordinate 
system in terms of which the polarization ratios of the two ellipses 
become complex conjugates of each other after a proper differential 
phase delay is introduced along the axes. The X’-Y’ coordinates cor- 
respond to the point C’, and the proper differential delay is represented 
by A on the Poincaré sphere as shown in Fig. 3a. Let us write down 
the expression for the polarization ratio in terms of the X’-Y’ 
coordinates! 


; e 2 tan ai 
Pt = mi + tan’ a;) — (1 — tan? a,) cos 28% peas la —tan?ai) sin =a 
(1 + tan? a:) + (1 — tan? a,) cos 26; (6) 
0 < |Pi < a when tana; > 0; 
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\ 
eae [| 


a las : lV CN 
CK, =2£, CK, =28', M{N=NM’, 
C CN 

(a) CK, =2B, C’K, = 2B'5 (b) 


Fig. 3—Simultaneous transformation of two polarization ellipses into two ellipses 
with parallel axes and equal axial ratios. 


where 7 = 1, 2; 


T T 
=e ge whem ein 2Be 270s 


Combining the condition |P{| = |P3| and the relation 63 = 6; — 8, 
one obtains 
(1 + tan? a2) (1 — tan? a1) 
Fe ee ee COS: 20 
1 (1 + tan? a,)(1 — tan? a2) 


Tv 
’ = — tan" —|, 0<~P <- (7 
, 2 sin 20 1 2 7) 





which fixes the X’-Y’ axes. The required differential phase delay A 
along the Y’ axis is determined by 


IP) — A = — (|P$— A) 
which gives 
A = 4([P{ + |PS). (8) 


The above differential phase shift corresponds to the rotation of the 
great circle arc My} iM, about the point © C’. The transformed d ellipses are 


represented by M; and M3 where 1 MIN = MIN and MiMi is per- 
pendicular to the equator as shown in Fig. 3b. The parallel major axes 
of the transformed ellipses are oriented with respect to the X’ axis at 
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TRANSMITTING | TRANSMISSION | RECEIVING 
ANTENNA ANTENNA 





(b) 


Ad,, dy DIFFERENTIAL PHASE — SHIFTER 
Aa,, Aa, DIFFERENTIAL ATTENUATOR 


Az PHASE — SHIFTER WILL PERFORM POLARIZATION TRACKING FOR ORTHOGONAL LINEAR 
POLARIZATIONS (A> PHASE — SHIFT AT BOTH TRANSMITTING AND RECEIVING ENDS WILL 
BE NEEDED FOR ORTHOGONAL CIRCULAR POLARIZATIONS) 


Fig. 4—Locations of the orthogonalization device. 


oN 

an angle equal to 4 NC’. Using eq. (13), we have 

2|Pi| 
1 — |P{|? 
Now two oppositely rotating circular polarizations can be obtained by 
imposing a differential attenuation along the major axis of the trans- 
formed ellipses. The value of differential attenuation is also given by 
tan (W/2) where y is determined by eq. (4) of the preceding section. 

In addition to the differential phase shift for orthogonalization as 
discussed above, circular polarization systems require A(7/2) phase 
shift* at both transmitting and receiving ends to convert to the final 
linearly polarized ports. The same amount of additional differential 
phase shift overall will also be needed in a system using orthogonal 
linear polarizations, if a Az phase shifter is used for polarization 
tracking. 


LN 
3 NC’ = $ tan” 


cos 3 [|Pi — |P2]}- (9) 


III. DISCUSSION 


The above analysis assumed that the device for orthogonalization 
would be located at the receiving terminal as shown in Fig. 4a. Since 
there always exist two elliptic polarizations which would become 
orthogonal after going through a linear transmission system with a 
certain polarization distortion, one can also put the differential ele- 


* The notation, A(z/2), implies a differential phase shift of +/2. 
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ments at the transmitting terminal as shown in Fig. 4b. Another ob- 
vious corollary states that the dual-polarization radiation of an an- 
tenna always can be orthogonalized in any particular direction. The 
differential attenuator should be located as illustrated in Figs. 4a and 
b at the receiver or the transmitter in order to satisfy the condition 
for minimum differential attenuation. Sometimes it is desirable to use 
differential elements at both the transmitting and receiving terminals. 
For example, one may wish to eliminate the polarization distortions 
of the transmitting and receiving antennas separately. 

It is often claimed that the use of circular polarization in a satellite 
communication system eliminates the need for polarization tracking. 
In order to realize this advantage, the depolarization of the satellite 
antenna radiation must be kept small over the entire coverage of ground 
stations. Furthermore, the matching requirement at each discon- 
tinuity of the waveguide feeding network and the radiating system is 
more stringent for circular polarization, because multiple reflections 
among the discontinuities often corrupt circular but not linear 
polarization. 
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APPENDIX 
Poincaré Spherical Representation 


A polarization ellipse shown in Fig. 5a is completely characterized 
by the axial ratio, A = minor axis/major axis, and the orientation of 
the ellipse. The sense of rotation can be taken into account by giving 
+ or — sign to the axial ratio. Now we define tan7! A as the ellipticity 
angle a with —45° S a S 45°, and take the angle between OX and 
the major axis as the orientation angle 6 with —90° < 6 < 90°. Then 
a point M on a sphere with longitude 28 and latitude 2a as shown in Fig. 
5b completely specifies a state of polarization. 

The points on the equator represent linear polarizations. The arc 
between two points along the equator is twice the angle between two 
linear polarizations. The points on the upper and lower hemispheres 
correspond to clockwise and counterclockwise (wave approaching) 
elliptical polarizations respectively, while the poles designate circular 
polarizations. For each point which represents an elliptic polarization, 
the projection K onto the equator defines the orientation of its major 
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(b) 


Fig. 5—Poincaré spherical representation of a polarization ellipse. 


axis. Any set of X-Y coordinates can be specified by its X direction 
which in turn corresponds to a point C on the equator. Then the 
longitude of KK can be measured with respect to C. 

The polarization ratio P is a complex number defined as the ratio 
between the Y and X components of the electric vector. Using the 
expression for the polarization ratio in eq. (6) and the following for- 
mulas for the spherical triangle shown in Fig. 5b, 


cos 2y cos 2a cos 28 (10) 
tang = tan 2a csc 2, (11) 


one can identify P = tan ye’*. The orientation and the ellipticity 
can be expressed in terms of P as follows: 


sin 2a = sin 2y sin ¢ (12) 
tan 28 = tan 2y cos ¢. (13) 


A differential phase delay of A of the Y component with respect to 
the X component will correspond to a clockwise rotation A of the arc 
oN 


CM around the point ‘‘C’”’ on the Poincaré sphere. 
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Theory of Roundoff Noise in Cascade 
Realizations of Finite lmpulse 
Response Digital Filters 


By D. 8. K. CHAN and L. R. RABINER 
(Manuscript received September 14, 1972) 


This paper presents a theoretical treatment of the roundoff noise problem 
for the special case of cascade realizations of Finite Impulse Response 
(FIR) digital filters.' Explicit relations for evaluating roundoff noise with 
the usual assumption of uncorrelated samples are presented. Useful scaling 
methods are stated and classified as to conditions when these methods are 
optimum. Important differences between use of these scaling procedures for 
Infinite Impulse Response (IIR) filters and FIR filters are pointed out. 
Finally, useful properties of linear phase FIR filters are discussed. 


I. INTRODUCTION 


In recent years, many techniques have been developed for the 
design of Finite Impulse Response (FIR) digital filters.2-° It is now 
possible to readily design filters with arbitrary frequency or time 
response characteristics using the windowing,?? frequency sampling,‘ 
or optimal design®-* methods. While both the windowing and fre- 
quency sampling techniques yield suboptimal filters, they are useful 
because of their simplicity and ease of design. The optimal design 
technique is of special importance because the filters it generates can 
be proved to be optimum in a certain sense,’ and because efficient 
algorithms exist for its implementation.’”:® 

As a result of these important developments, the FIR type of 
digital filter is becoming increasingly attractive as an alternative to the 
IIR (Infinite Impulse Response) type of filter for practical applications. 
A major advantage of FIR filters over IIR filters is that an FIR filter 


t This paper is based on a thesis! submitted in partial fulfillment of the require- 
ments for the degrees of Bachelor of Science and Master of Science in the Department 
of Electrical Engineering at the Massachusetts Institute of Technology in September 
1972. 
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can have an exactly linear phase response while approximating an 
arbitrary magnitude frequency response. But even without considering 
this important advantage, current research’ is revealing that in certain 
cases FIR filters are competitive with IIR filters in terms of speed and 
cost. Thus the implementation of FIR filters using finite-precision 
arithmetic is becoming an important area for research. 

Up to the present, little is known as to how different types of FIR 
filter realizations behave with respect to quantization effects. Since 
hardware, specifically for the purpose of realizing FIR filters, has 
already been built for experimentation by various research groups,!?"!! 
it is important to obtain more knowledge to guide the implementation 
phase of FIR digital filter design. The purpose of this paper is to 
present a theoretical treatment of several problems associated with 
implementations of these filters. 


II. PRELIMINARY REMARKS 


The effects that quantization has on an IIR filter can be classified 
into three basic categories: 


(t) Quantization of the values of samples derived from a con- 
tinuous input waveform causes inaccuracies in the representa- 
tion of the waveform (A-D noise). 

(ii) Finite-precision representation of the infinite-precision filter 
coefficients alters the frequency response characteristics of the 
filter (coefficient accuracy problem). 

(iz) Finite-precision arithmetic causes inaccuracies in the filter 
output (roundoff noise) which, together with the finite dynamic 
range of the filter, limit the signal-to-noise ratio attainable. 
Also, finite-precision arithmetic can lead to limit cycles where 
the output samples are generally highly correlated. 


These same quantization effects also occur in finite wordlength FIR 
filters with the important exception that limit cycles cannot occur in 
nonrecursive realizations of FIR filters. In this paper only the third 
type of quantization effect, viz., roundoff noise, will be discussed. 
Furthermore, fixed-point arithmetic with rounding will be assumed. 

Except for the first category above (A-D noise), all quantization 
effects depend in degree and character on the type of structure used to 
implement a filter. There are three well-known structures in which an 
FIR transfer function can be realized. They are the direct form, the 
cascade form, and the frequency-sampling structure.’? Other less 
well-known structures based on polynomial interpolation formulas 
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= Vm Vr 
o © 


Fig. 1—Cascade-form filter section. 


are also possible; these include the Lagrange, Newton, Hermite, and 
Taylor structures. However, it is as yet unclear under what circum- 
stances, if any, these other structures may be more advantageous than 
the well-known structures. 

Only the cascade structure will be discussed in this paper. A second- 
order filter section, as shown in Fig. 1, will be used as the basic building 
block for the cascade structure. Although several minor variations to 
this configuration for the filter sections are possible,! the results pre- 
sented here are sufficiently general so that they can be readily applied 
to other configurations as well. 

Aside from section configuration, the prime issues that must be con- 
fronted in the realization of filters in cascade form are scaling and 
section ordering. Proper scaling must be performed on a filter in 
cascade form in order that full use of the dynamic range of each section 
can be made while avoiding error-producing overflows. By proper 
scaling, the signal-to-noise ratio of a filter can be maximized for a 
given quantization step size and section ordering. 

Proper ordering must also be determined for a filter in cascade form 
if the filter is to be useful at all, since the noise output of a cascade 
filter can depend dramatically on the way it is ordered. For example, 
Schtissler® showed a 32nd-order FIR filter which, ordered two different 
ways, yields output noise variances that differ by a ratio on the order 
of 10°. The problem of section ordering for cascade FIR filters has been 
investigated in depth,!' and the results show that for higher-order 
filters, the variation of output noise variance across all orderings is 
much greater than 10°. 

Jackson!> has formulated the roundoff noise problem for a general 
digital filter and has proposed an approach to the scaling of filters to 
satisfy dynamic range constraints. Most of his results can be specialized 
to the case of FIR filters by assuming a constant polynomial in the 
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denominator of the transfer function. However, the different perspec- 
tive obtained by studying the FIR case separately affords many 
additional insights. 

In this paper, formulas for the evaluation of roundoff noise variances 
in the FIR cascade structure are presented. Also, specific scaling 
methods for FIR filters are defined, two of which can be proved to be 
optimum for two classes of input signals. Finally, certain properties 
of the linear phase cascade structure useful in the study of section 
ordering are stated and can be proved. However, for reasons of space, 
no proofs are included in this paper. They can be found in Ref. 1. 


III. DEFINITIONS 


The general transfer function for an N-point FIR filter can be 
written in the form 


N-1 
H(z) = Xo h(k)z* (1) 
k=0 
where the real-valued sequence {h(k), k = 0, ---, N — 1} is the im- 


pulse response of the filter. Alternatively, H(z) can be expressed in the 


factored form 
Ns 
H(z) = TI (bo: + bie + bei2~*) (2) 
t=1 
where b;;, 7 = 0, 1, 2,7 = 1, ---, Ns are real numbers and N,, the 
number of factors, is defined as 


N-1 
N odd 





N, = 
N 
— N even 
2 
and bey, = 0 if N is even. 
A linear phase filter is defined to be a filter, the transfer function H (z) 
of which is expressible in the form 


H (2) | sme” = Hei) = +|H (e%) |e7s0# (3) 


where a is a real positive constant with the physical significance of 
delay in number of samples. The factor + is necessary since H (e”) 
actually is of the form 


H(e%) = H*(e)e—see 
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where H*(e%“) is a real function taking on both positive and negative 
values. It is useful to define a mirror-image polynomial (MIP) of 
degree N to be a polynomial of the form ~~ a2", the coefficients of 
which satisfy the relation 


Qp = Aan—k OSKEN. 


Necessary and sufficient conditions on H(z) such that a filter with 
transfer function H(z) has an exactly linear phase response can then 
be stated as follows: 


Theorem 1: H(z) can be expressed in the form (8) if and only if one of 
the following equivalent conditions hold: 


(1) h(kh) =hA(N-—-1-—k), 085k 5 N—1. 
(ii) If 2; is a zero of H(z), then 2; ' is also a zero of H(z). Also, if 
2: = +178 a zero of H(z), then it occurs in even multiplicity. 
(tit) Suppose 2; ts a zero of the ith factor in (2). Let S = {1: 2; ts real} 
andQ = {1:7 E S}. Thenf(z) = Ties (boi + biz! + boiz?) ts 
a mirror-image polynomial in 2, and for all 1€CQ, either 
bo: = be; or there exists 7 ¥ 1, 7 C Q, such that 


ea : 


Furthermore, the following is a sufficient condition for H(z) to be expres- 
sible in the form (8): 


(wv) In (2), for1 SiS Ng, either bo; = Oand bo; = b1:, or boi = dai, 
or there exists 7 ¥ 1,1 5 9 S Na, such that 


In all cases the value of a isa = (N — 1)/2. 


It should be pointed out that a section with bo; = be; is necessarily 
one which synthesizes either two complex conjugate zeros on the 
unit circle, or two reciprocal zeros on the real axis, or two identical 
zeros at +1 or —1. Furthermore, two sections which satisfy (4) are 
precisely those sections which synthesize reciprocal zeros (i.e., if 2; 
is a zero of one section, then z; ‘ isa zero of the other section). Thus, 
taking (2) as the basis for the FIR cascade form, condition (iv) of 
Theorem 1 provides a way to assign zeros to individual sections of 
the cascade structure so that linear phase is guaranteed independent 
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of scaling or ordering. Hence, in this paper the following convention 
of zeros assignment for linear phase filters will be adopted: complex 
zeros are grouped by conjugate pairs, real zeros that are reciprocals of 
each other are paired together, while doubled or higher multiplicity 
zeros are grouped by pairs of the same kind. In this way the only zero 
that can occur by itself in a section is z = — 1 (since by Theorem 1, 
z = + 1 is not allowed as a zero of odd multiplicity). 

The definition (8) of a linear phase filter requires the filter to have 
both constant group delay and constant phase delay. However, if 
only constant group delay is desired, a second type of ‘‘linear phase’’ 
filter can be defined in which the phase of H (e%) is a piecewise linear 
function of w, i.e., 


H(e%) = +|H (e%) | e96-ae), (5) 


It can be shown! that with the constraint (1) on the form of H(z), the 
only possible solutions for 8 e[—7z, 7] is 8B = + (kr/2), k = 0, 1, 2. 
If 8 = 0, +7, (5) reduces to (8). Thus the only new cases added are 
when 6 = + 7/2. These cases arise exactly when z; = + 1 occurs as 
a zero of H(z) in odd multiplicity, or, equivalently, when {h(k)} 
satisfies 


h(k) = —h(N-1—-hk) OSKkEN-1. 


Filters of this special type are useful in the design of wideband differ- 
entiators.!° However, this type of filter will not be considered in this 
paper and the term “‘linear phase filter’’ will be restricted to refer to 
those filters satisfying (3). 


IV. THEORY OF FIR CASCADE STRUCTURES 
4.1 Roundoff Noise in the Cascade Structure 


The analysis of roundoff noise in this paper is based on the usual 
model used for such analyses in digital filters.15:1718 In particular, each 
multiplier in a filter is modeled as an infinite-precision multiplier 
followed by a summation node where roundoff noise is added to the 
product so that the overall result equals some quantized level. Each 
noise sample is modeled as a random variable with uniform prob- 
ability density on the interval (—Q/2, Q/2) and zero density elsewhere, 
where Q is the quantization step size. Thus each sample is a zero-mean 
random variable with a variance of Q?/12. 

Furthermore, the following assumptions are made: 


(1) Any two different samples from the same noise source are 
uncorrelated. 
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(11) Any two different noise sources (i.e., associated with different 
multipliers), regarded as random processes, are uncorrelated. 
(zit) Each noise source is uncorrelated with the input signal. 


Thus each noise source is modeled as a discrete stationary white 
random process with a uniform power density spectrum of magnitude 
Q?/12. 

Applying this model to the filter section shown in Fig. 1, the addition 
of a noise source to the output of any multiplier is seen to be equiv- 
alent to adding a noise source to the output of the section. Therefore, 
to model a section of a cascade filter, k; noise sources are added to the 
output of the section, where k; is the number of multipliers with non- 
integer coefficients in the section. Or, equivalently, by assumption (72) 
above, one noise source of variance k;(Q?/12) can be added instead. 

For the configuration shown in Fig. 1, k; is in general equal to 3. 
However, when bo; = b2;, the signals of the two branches feeding the 
multipliers with coefficients bo; and b:; can first be summed before 
being multiplied by the common coefficient, thus reducing k; to 2. 
Furthermore, by a sacrifice in speed (assuming serial arithmetic), 
it is possible, as demonstrated by practical hardware,” to reduce 
k; to 1 for all 7 by summing all products in each section before per- 
forming rounding. It is of interest to point out that the same can 
be done in the direct form, resulting in effectively only one noise 
source of variance Q?/12 feeding into the output of the filter. 

Before proceeding further, some notations need to be developed. 
Let H;(z) denote the transfer function of the ith section of a filter 
H(z), ie., 


Ns 
A(z) = II H(z) (6) 
i=1 
where 
A j(z) = bo: + biz + dee. (7) 


As a convention, filter sections will be numbered in increasing numbers 
according to increasing distance from the filter input (i.e., the section 
at the input is called the Ist section). 

Furthermore, define 


Ns 
Il 4;(z) 0OsS718SN,-1 
Gi(z) = 4 jor41 : (8) 
1 i= N; 
and let {g:(k)} be the impulse response of G;(z), i.e., 


Gilz) = Le gi(k)e*. (9) 
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e, (n) @9 (n) CN, (n) 





Fig. 2—Equivalent models for a filter in cascade form. 


Then a cascade filter can be modeled as in Fig. 2a or equivalently 
as in Fig. 2b, where {e;(n)} is the noise source for the 7th section. 
Letting {f£;(n)} denote the noise sequence at the filter output due to 
the zth noise source alone gives 


E,(n) = DV gilk)e(n — k). (10) 


By the stationarity of {e;(n)}, the variance of E;(n) is independent 
of n; hence denoting this variance by o7, assumption (7) above leads 
to the relation 


oi = Ngilke(n — ky? 


Q?_ 2 
io (k) ) 
Now the total noise output is given by 


E(n) = & En) = EE gehen — b). (12) 


t=1 k 
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Therefore by assumptions (7) and (72) 


pes) Bs a. (13) 


t=1 


4.2 Methods of Scaling to Meet Dynamic Range Constraints 


A practical digital filter, necessarily implemented as a physical 
device, must have a finite dynamic range. Especially when fixed-point 
arithmetic is employed, this dynamic range sets a practical limit to 
the maximum range of signal levels representable in a filter and acts 
to constrain the signal-to-noise ratio attainable. 

In some filter structures such as the direct form, given the filter 
transfer function, the designer has no control over the relative signal 
levels at points within the filter. Only the gain of the overall filter 
can be varied. However, in a cascade realization with N, sections 
there are N, — 1 degrees of freedom available in addition to the 
overall filter gain and the ordering of sections. 

To see this, a factorization for H(z) is defined which is unique up 
to ordering of factors, in the form 


H(e) = 6 IL Aa) 


A(z) = ao + aye + are? (14) 


where {a,;} satisfies 


M x 


Aoi = 0, 


lax] =1 t¢=1,---, Ny. (15) 


0 


I 


2 


Then the transfer function for the 7th section in a cascade realization 
can be written as 


H;(z) => S:H:(z) (16) 

where S; is an arbitrary constant, subject only to the constraint that 
Ns 

IIS; = 8. (17) 
i=1 


Thus given 6, VN, — 1 of the S,’s can be chosen at will. 

Any rule for assigning values to {S;} will be referred to as a scaling 
method. Obviously, some scaling method must be employed in the 
design of a cascade filter whether or not one is concerned with dynamic 
range constraints, since numerical values must be assigned to the S;’s. 
When dynamic range is an issue, the constraints it imposes can be 
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met in some best manner by choosing the proper scaling method. 
In this paper, filters, designed so that no arithmetic overflow in them 
can cause distortion in the filter output, will be studied. Therefore, the 
investigation of scaling methods will be restricted to those methods 
which guarantee that for a given class of input signals no distortion- 
causing overflow occurs in the scaled filter. 

It can be shown" that, in an addition operation, if two’s comple- 
‘ment arithmetic is used, as is usually the case, then as long as the final 
result is within the representable numerical range, individual partial 
sums can be allowed to overflow without causing inaccuracies in the 
result. In this paper, it is assumed that all additions in a filter are 
done using two’s complement arithmetic. Then, to guarantee that no 
distortion caused by overflow occurs at a cascade filter’s output, only 
the input and output of each filter section need be constrained not to 
overflow. 

To simplify the discussion of scaling methods, the following defini- 
tions are used. Let 


Fila) = & flke* = TT Hye) (18) 
and = a P2eeN; 
Fiz) = X flkye* = IT Ai). (19) 


Also, let {v;(n)} be the output sequence of F(z) or H;(z). Furthermore, 
assume that the maximum magnitude of numerical data representable 
in a filter is 1.0. Then the necessary overflow constraints on a cascade 
filter can be stated as 


lu(m)|S1 1<iSQN,, alin. (20) 


Necessary and sufficient conditions for (20) to hold for two classes 
of input signals are given below. Theorem 2 deals with the class of 
input sequences {x(n)} which satisfies |x(n)|S 1 for all n. For sim- 
plicity, this class is refererd to as “‘class 1.’”’ Theorem 3 deals with 
the class of inputs with transform X (e%”) which satisfies 


1 Qn 
~| | X(e%) |dw <1. 
Qn 0 
This class will be called ‘“‘class 2.”’ By virtue of the fact that 


r= a X (e7*)eF¢"™dw (21) 
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and hence 
1 27 
Ix(n)| <— if |X (ei) |de, (22) 
27 0 
class 2 is a subset of class 1. 


Theorem 2: Suppose |x(n)|S 1. Then |v:i(n)|S 1, 1 SiS N, and 
all n, tf and only if 


IA 
— 
) 
I 
tS 
= 
&% 


2 
DV |fth)| S (23) 
k=O 

Theorem 8: Suppose 1/(2r) Jo*™ |X (e%)|dw <1. Then |v;(n)| S 1, 
1sisN, and all n, tf and only if 


IFi(e*)|S1 G¢=1,---,N, OSwS 2. (24) 
Conditions (23) and (24) of Theorems 2 and 3 can be restated to 
give conditions on {S;}. Recall that the H;(z)’s are unique once H (z) 


is given, hence the F;(z)’s and {f;(k)}’s are also unique. Equations 
(16), (18), and (19) give 


fuk) = (IL Ss) fle (25) 
and * 
ria) = (IL Si) A). (26) 
rm 
Therefore, conditions (23) and (24) can be restated respectively as 
t 20 —1 
isi s[ = 140i | (27) 
and : e a=1,--::,N,. 
IT |Si| $C max | F(e%) | > (28) 
I= =o =27 : 


These then are conditions which, for the class of inputs concerned, a 
scaling method must satisfy. It will be shown next that in some sense 
optimum scaling methods are obtained when (27) and (28) are satis- 
fied with equality. For ease of reference, two scaling methods will 
first be defined. 

Define sum scaling to be the rule: 


Ws=[5 imi] i=, (29) 
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or stated recursively, 
PCa 
8: = a 2% -1 (30) 
[(1s:) = Awl fo b= we 


Also, define peak scaling to be the rule: 


ID S;=[ max {A(e*)|}2  t=1,---,N, (31) 
j=1 0 Sw S2r 
or a 
CL max |Fi(e%)|[JT'  i=1 
Se 2 (32) 
| (11 s;) max |Flc)| | i= 2,---,Ns. 
j=1 0Sw S2r 


Theorem 4: Given an FIR transfer function to be realized in cascade 
form (as defined in Fig. 1) using fixed-point arithmetic of a given word- 
length, and given the ordering of filter sections, assume that: 


(t) The number of noise sources in each section (t.e., ki) is inde- 
pendent of the scaling method. 
(it) All filter coefficients can be represented to arbitrary precision. 
(i2t) No overflow is allowed to occur at the input and output of each 
section. 
(zv) The overall gain of the filter is maximized subject to no overflow 
at the filter output. 


Then each of the following scaling methods is optimum for the class of 
input signals stated, in the sense that it yields the minimum possible 
roundoff notse variance as defined in (18) among all scaling methods 
which satisfy conditions (211) and (iv) above for the class of inputs 
considered. 


(t) Sum scaling for class 1 signals. 
(iz) Peak scaling for class 2 signals. 


Thus optimal scaling methods are established for two classes of 
input signals. It is possible to define other classes of signals by con- 
sidering the “Lp norm’’ of their transforms.'®” Specifically, the Lp 
norm of X (e’) is defined as 


1 Qr 1l/p 
|X(e*)]|» = | = | |X(o) [rao | L<p<~ (33) 
To 
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where for p = © the limit as p— © of the right-hand side is meant. 
For each p, a class of signals can be defined consisting of those se- 
quences with transforms which satisfy 


|X (e%) ||» 1. (34) 


Signals satisfying (84) will be referred to as Lp-norm constrained 
signals. Note that Zi-norm constrained signals are simply class 2 
signals. 

For proofs of the following useful theorem, refer to Refs. 1, 20, and 21. 


Theorem 5: Let X (e%) and Y (e%) be transforms of sequences. Then 


(2) |X(e%) ||. = max |X(e%)| 
(it) |X (e%) Y(e%) ||, S ||X(e*)|| oll ¥(e*) lla 
uf 1/p+1/q=1, lspqso 
(iit) |X (e%) ||, S |X(e*) ||, if lSrS<sso. 


Since with input {z(n)}, 


= 1 jo) X (es) esjond 
os(n) = = : (ei) X (e##)e30%des, (35) 
so that : , 
lon)| Ss | [FP (e%*)X(e) [dw = ||F.(e%*)X(e%)|la, (36) 


by Theorem 5 (22), 
|vi(n)| S||Fs(e**) || |X (e*)|[¢ 1 


IIA 


tS N.. (37) 


Hence for Z,-norm constrained signals, i.e., if ||X(e*)||, S$ 1, the 
following scaling method (Z,-norm scaling) is obtained. 


q 
Pp =>—=_l 
[Fe*)|lp = 1 = (38) 
es i=1,--',N,, 
or stated in terms of {S;}, 
ILS; = [l|F(e*)pf* t= 1,---, Ns. (39) 
j=l 


Notice that by virtue of part (¢) of Theorem 5, Z,.-norm scaling is 
just peak scaling which has been shown to be optimum for class 2, or 
Z,-norm constrained, signals. Furthermore, by part (777) of the theorem, 
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the following hierarchy of classes of signals is obtained : 


class 1) class 2 L,-norm constrained D L,-norm 
constrained 
filspsqs~. 


In general, class 1 and class 2 signals are the most useful to consider. 
L,-norm constrained signals with Le-norm scaling are useful when all 
inputs to a filter have finite energy bounded by a known value. For, by 
Parseval’s Theorem, the energy of {x(n)} is simply given by 
(|| X (e4)||2)?. Hence, if the input signals are first scaled so that their 
maximum energy is 1.0 (or the squared dynamic range of the filter), 
then L».-norm scaling is sufficient to ensure no overflow. 

L.-norm scaling finds greater application for FIR filters than for 
IIR filters because, in the former case, it is applicable for a larger class 
of input signals. In particular, an Nth-order FIR filter has only N 
samples of memory; thus if the input signal to an Nth-order IIR 
filter consists of bursts of energy spaced more than N samples apart 
with zero energy in between, then the filter will effectively ‘‘see’’ only 
one burst at a time. Hence, the maximum energy of a burst can be 
used as the bound on the energy of the input as far as scaling is con- 
cerned. Thus an infinite-energy signal can have the effect of a finite- 
energy signal on an FIR filter. 

Clearly, sum scaling and peak scaling can also be applied to IIR 
filters.'5 In fact, Theorems 2 and 3 are also valid for IIR filters. 
However, the input sequence needed in Theorem 2 to prove necessity 
in the case of IIR filters is an infinite-duration sequence extending to 
— © with full dynamic range magnitudes, and signs that match those 
of { f:(k)} for some 7. Since { f;(k)} for IIR filters is infinite in duration 
for all 7, clearly such an input sequence is highly improbable. Hence, 
class 1 signals have been deemed too restrictive a description for 
ordinary inputs to an IIR filter, resulting in too stringent a scaling 
method.!5 

However, for FIR filters it is not difficult to find an input sequence 
within dynamic range which will require sum scaling to ensure no 
overflow, since only a small, finite portion of the sequence need match 
up with the {f;(k)}’s. For example, if F:(z) has a zero with angle wo, 
a/2 S wo <7, then all three samples of { f1(k)} have the same sign; 
hence an input sequence need only have three consecutive samples of 
value 1 before |v:(m)| = >2.|f1(k)| for some n. 
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4.3 Properties of the Linear Phase Cascade Structure 


Two theorems regarding certain properties of the linear phase 
cascade form are now given. These results are useful in the investiga- 
tion of ordering of cascade filter sections." 


Theorem 6: Let H;(z) be the transfer function for the ith section of a 
linear phase FIR filter in cascade form, where 


H(z) = bo: + bie + boiz, 
and let w; be the angle of one of its zeros, —7 Sw; S a. Then for allt: 


|H(e*)| 0S Jail < 5 
(2) max |H,(e)| = 


IA 
FS) 


| H :(e) | < |o;| 


2 
(72) = | bis| = max | H:(e%*) | = max (| H:(e*°) ; | H,(e?7) | ) ‘ 


The next theorem is concerned with the equivalence of certain order- 
ings with regard to output noise variance. In particular, it states that 
with peak scaling each pair of sections in a filter which synthesize 
reciprocal zeros is completely interchangeable without affecting the 
output noise variance of the filter. With sum scaling, however, this 
is not necessarily true. Nevertheless, a weaker condition can be stated 
which says that, with sum scaling, if every pair of sections which 
synthesize reciprocal zeros of a filter is interchanged in position, then 
output noise variance is not changed. Figure 3 illustrates two such 
equivalent orderings. 


Theorem 7: Let {H;(z)} and {H;(z)} be the section transfer functions of 
two orderings for a linear phase filter H (z), both scaled by the same method, 
thus 


Ns Ns 
A(z) = JI Hi(z) = I] Ailz). 
t=1 i=1 
Then filters with section transfer functions {H;(z)} and {H;(z)} produce 


identical output noise variances if either of the following conditions ts 
true: 
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40 


O6 


Re [z] 


O6 





Fig. 3—Two orderings with equal output noise variances. 








(i) Peak scaling is used and for each i, H;(z) and H;(z) have either 
the same zeros or reciprocal zeros [i.e., H;(z:) = Oif H;(zz') = 0). 

(ii) Sum scaling is used and for all 7, 2;' is a zero of H;(z) whenever 
2; 1s a zero of H(z). 
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An Algorithm for Minimizing Roundoff 
Noise in Cascade Realizations of Finite 
Impulse Response Digital Filters 


By D. 8. K. CHAN and L. R. RABINER 
(Manuscript received September 14, 1972) 


Experimental results on roundoff noise in cascade realizations of 
Finite Impulse Response (FIR) digital filters are presented in this paper.* 
The entire roundoff noise distribution (i.e., over all possible orderings) 
is given for several low-order filters using both sum and peak scaling. 
Based on observations about this distribution, as well as intuitive argu- 
ments about the effects of ordering on roundoff noise, an algorithm for 
minimizing roundoff noise 1s presented. Experimental verification of this 
algorithm for a wide range of filters 1s given. 


I. INTRODUCTION 


As discussed in previous works,!* the implementation of FIR filters 
using finite precision arithmetic has become an important issue in 
recent years. For cascade realizations of FIR filters, roundoff noise 
is a crucial problem. In Refs. 1 and 2, some of the theoretical bases for 
the analysis of roundoff noise in the FIR cascade form have been 
considered. This paper presents a large body of experimental results 
which depict the dependence of roundoff noise on several of the 
important parameters of a cascade FIR low-pass filter. Most impor- 
tantly, these results point to an algorithm which can find efficiently, 
for a cascade filter, an ordering which has a noise variance very close 
to the minimum possible. Experimental verification of this algorithm 
for a wide range of filters is presented. 

Low-pass, extraripple? filters are used throughout these investigations 
as being representative of FIR filters. It will be seen that most results 

* This paper is based on a thesis! submitted in partial fulfillment of the require- 
ments for the degrees of Bachelor of Science and Master of Science in the Department 


of Electrical Engineering at the Massachusetts Institute of Technology in September 
1972. 
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should not depend on the type of filter used. Figure 1 shows the mag- 
nitude response of a typical extraripple filter (which by definition has 
linear phase) and the parameters which define it. Of the four parameters 
F,, F2, Dy, and Ds, only three can be independently specified. The 
parameters N (impulse response length), Fi, D1, and De» will be chosen 
as independent variables in these investigations. The studied ranges of 
variation of these parameters are asfollows:7 < N S 129,0.12 D,2 
0.001, 0.1 = D. 2 0.001, 0 < F, < 0.5, and 0 < F, < 0.5 (sampling 
frequency = 1). These ranges comprise a large range of the signifi- 
cant values that these parameters take on. In the present state of the 
art in real-time digital filter hardware, 128th-order (NV = 129) is the 
highest order that can be implemented in cascade form at a sampling 
rate of 10 kHz (e.g., typical for speech processing) .*»> Furthermore, the 
stated ranges for D,; and Dy, are those significant for many speech 
processing systems.°® 

While a great deal of experimental data has been collected, only 
representative examples will be presented here. For more examples 
see Ref. 1. 


II. PRELIMINARY REMARKS 


In Ref. 2 it is shown that given a transfer function H (z) to be realized 
in cascade form and the order in which the factors of H(z) are to be 
synthesized, there remain N, degrees of freedom (including gain of 
filter) in the choice of filter coefficients, where N, is the number of 
sections of the filter. Scaling methods are developed to fix these N, 
degrees of freedom, and two particular methods, viz., sum scaling and 
peak scaling,* are shown to be optimum for the particular classes of 
input signals which they assume. These scaling methods will be 
applied in this paper. 

The prime issues in the realization of filters in cascade form are 
threefold-scaling, ordering, and section configuration. Because of the 
simplicity of a 2nd-order FIR filter, there is little freedom in the choice 
of a structure for the sections of a cascade filter. In Ref. 2, the con- 
figuration shown in Fig. 2a is assumed because it turns out to be the 
most useful in a practical situation. Another configuration, Fig. 2b, 
is also mentioned in Ref. 2. Because, as seen from Fig. 2c, these two 
configurations can be readily accommodated in a more general sub- 


*Sum (peak) scaling is defined to be a method of scaling where the scale factors 
for the cascade sections are chosen so that the sum of the impulse response magnitudes 
(peak of the magnitude of the frequency response) up to that section does not exceed 
one. 
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N=11 


F, = 0.216 


Fy = 0.284 
D, =0.1 


Dy =0.1 


MAGNITUDE 





0 0.1 0.2F, H Fo 0.3 0.4 0.5 
FREQUENCY 


Fig. 1—Definition of filter parameters D1, De, Fi, and Fo. 


filter structure, it is here assumed that the configurations of Figs. 2a 
and b are both used in a cascade structure, depending on whether 
boi % be: or bo; = be; respectively. The configuration of Fig. 2b has 
the advantage of having lower noise than the configuration of Fig. 2a. 
The option of summing all products in an increased length register 
before rounding is also possible for all configurations. However, the 
gain in signal-to-noise ratio does not seem to be worth the required 
sacrifice in speed (assuming serial arithmetic).”? In any ease, all re- 
sulting noise variances would simply be scaled down by a factor of 
from 2 to 3 if this strategy were used instead of rounding after each 
multiplication, as assumed here. 

Other possible section configurations will be discussed later on. Since 
scaling is treated in depth in Ref. 2, the major concern here is the 
ordering of sections. Unlike the scaling problem, no workable optimal 
solution (in terms of feasibility) to the ordering problem has yet been 
found for cascade structures in general. The dependence of output 
roundoff noise variance on section ordering, given a scaling method, 
is so complex that no simple indicators are provided to assist in any 
systematic search for an ordering with lowest noise. Any attempt to 
find the noise variances for all possible orderings of a filter involves on 
the order of NV! evaluations, which clearly becomes prohibitive even 
for moderately large values of N;. Thus there is little doubt that 
optimal ordering with time constraint is by far the most difficult 
issue to deal with in the design of filters in cascade form. 
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Fig. 2—(a) Cascade form filter section. (b) Alternate cascade filter section. (c) 
General cascade filter section. 


Since finding an optimal solution to the ordering problem through 
exhaustive searching is very time consuming (if not impossible in any 
feasible amount of time) for all but very low-order filters, it is im- 
portant to find out how closely a suboptimal solution can approach 
the optimum and how difficult it would be to find a satisfactory sub- 
optimal solution. Even this concern, however, would be unfounded if 
the roundoff noise level produced by a filter were rather insensitive to 
ordering. Then the difference in performance between any two order- 
ings may not be sufficient to cause any concern. However, Schiissler 
has demonstrated that quite the contrary is true.”’ He showed a 33- 
point FIR filter which, ordered one way, produces o? = 2.4 Q? (where 
Q is the quantization step size of the filter), while ordered another way 
ylelds o? = 1.5 X 108 Q? (assuming all products in each section are 
summed before rounding). This represents a difference of 1.6 bits versus 
14.6 bits of noise. Clearly, the difference is large enough so that the 
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problem of finding a proper ordering of sections in the design of a cas- 
cade filter cannot be evaded. 

An important question to pursue in investigating suboptimal solu- 
tions is whether or not there exists some general pattern in which 
values of noise variances distribute themselves over different orderings. 
For example, for the 33-point filter mentioned above, are all noise 
values between the two extremes demonstrated equally likely to occur 
in terms of occurring in the same number of orderings? Or, perhaps, 
only a few pathological orderings have noise variance as high as that 
indicated. On the other hand, perhaps only very few orderings have 
noise variances close to the low value, in which case an optimum solu- 
tion would be very valuable, while a satisfactory suboptimal solution 
may be just as difficult to obtain as the optimum. 

In the next section, these questions will be answered by investi- 
gating filters of sufficiently low order so that calculating noise variances 
of N,! different orderings is not an unfeasible task. The implications of 
results obtained will then be generalized. 


III. CALCULATION OF NOISE DISTRIBUTIONS 
3.1 Methods 


The definitions of sum scaling and peak scaling in Ref. 2 indicate 
that, for FIR filters, sum scaling is much simpler to perform than peak 
scaling. To achieve peak scaling, the maxima of the functions F;(e*)* 
must be found for all 7 given an ordering. Even using the FFT, this 
represents considerably more calculations than finding >>?"_,| fi(k) | 
for all 2 as is necessary for sum scaling. In the 33-point filter mentioned 
above, Schiissler used peak scaling on both the orderings. It will be 
shown in Section IV that, given a filter, peak and sum scaling yield 
noise variances that are not very different (within the same order of 
magnitude), and, in fact, experimental results indicate that they are 
essentially in a constant ratio to one another independent of ordering 
of sections. Hence the general characteristics of the distribution of 
roundoff noise with respect to orderings should be quite independent 
of the type of scaling performed. In order to save computation time, 
sum scaling will be used in these investigations. 

Returning to the question of section configuration, for Infinite 
Impulse Response (IIR) filters Jackson? has introduced the concept of 
transpose configurations to obtain alternate structures for filter sec- 


* F,(ei*) and {fi(k)} are defined in Ref. 2 and are the frequency response and 
impulse response of the cascade of sections 1 to 7. 
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3—(a) Transpose configuration of Fig. 2a. (b) Alternate configuration of 
2b. 


Fig 
Fig. oa (c) Alternate configuration of Fig. 


tions. However, the application of this concept to Fig. 2a yields the 
structure shown in Fig. 8a, which is seen to have the same noise 
characteristics as the structure in Fig. 2a since, by the whiteness 
assumption on the noise sources, delays have no effect on them. There- 
fore, the structure of Fig. 3a has no advantages over other structures 
as far as roundoff noise is concerned. The only other significant alter- 
nate configuration for Fig. 2a is shown in Fig. 3b. The counterpart 
for Fig. 2b is Fig. 3c and is valid when bo; = be;. Both of these new 
configurations have exactly the same number of multipliers as the 
original ones. However, one noise source is moved from the output 
to essentially the input of the section. Thus it is advantageous to use 
the structures in Figs. 3b and c for the 7th section when 


1 
pr ILC) < E oith) (1) 


where {g:(k)} is, as defined in Ref. 2, the impulse response of the 
equivalent filter seen by the 7th noise source. However, in order to 
have no error-causing internal overflow when the input and output 
of a section are properly constrained, the structures of Figs. 3b and ¢ 
can be used only when by; S$ 1. If bo; > 1, either four multipliers are 
required, or Fig. 3b reduces to Fig. 2a. 

In the investigations that follow, for each section of a filter, the 
configuration among Figs. 2a, 2b, 3b, and 3c which is applicable and 
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Fig. 4—-Flow chart of scaling and noise calculation subroutine. 
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results in the least noise will be employed. It turns out that this 
flexibility in the choice of configuration has little effect on the noise 
distribution characteristics of a filter. For low-noise orderings, the con- 
figurations of Figs. 2a and b are almost always more advantageous 
than the other configurations. For high-noise orderings, the alternate 
configurations help to reduce the noise variance, but the difference is 
comparatively small. Thus in actual filter implementations the 
structures in Figs. 3b and c may be ignored. 

Figure 4 is the flow diagram of a computer subroutine which is used 
to accomplish scaling, choice of configuration, and output noise vari- 
ance calculation given a filter and its ordering. The input to the sub- 
routine consists of N, (the number of sections) and the sequence 
{C;, 1S j 54, 15725 N,}, the elements of which are unscaled 
coefficients of the filter, defined by 


A(z) = Cis(Coi + Cae + Caz) 1s71 


IIA 


N; (2) 


where H;(z) is the 7th section in the filter cascade. The sequences 
{fi(k)} and {g;(k)} in Fig. 4 are the impulse responses of the cascade 
of the first 7 sections and the last (NV, — 7) sections respectively. The 
coefficients {C;;} on input are assumed to be normalized so that, for 
all z, Ci; = 1 and at least one of Co; and C4; equals 1. On return {C;;} 
contains the scaled coefficients and Nx is the value of output noise 
variance computed in units of Q?, where Q is the quantization step 
size of the filter. 

Using this subroutine, the noise output of all possible orderings of 
several FIR filters ranging from N, = 3 to N, = 7 was investigated. 
By Theorem 7 (77) in Ref. 2, for any filter with at least one set of two 
complex conjugate pairs of reciprocal zeros, there are at most N,!/2 
orderings that differ in output noise variance. This is true since if all 
orderings are divided into two groups, according to the order in which 
the reciprocal zeros are synthesized in the cascade, then Theorem 
7(iz) establishes a one-to-one correspondence between each ordering 
of one group and some ordering of the other group. Thus, in the investi- 
gation of all possible noise outputs of a filter, where possible, a pair 
of sections which synthesize reciprocal zeros is chosen and all order- 
ings in which a particular one of these sections precedes the other are 
ignored. 


3.2 Discussion of Results 


Using the methods and procedures described, the noise distributions 
of 27 different linear phase, low-pass extraripple filters were investi- 
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Fig. 5—Positions of the zeros of one filter. 


gated. Twenty-two of these filters are 13-point filters, since N = 13 
represents a reasonable filter length to work with. Thirteen-point 
filters have six sections each, corresponding to 6! or 720 possible order- 
ings of sections. By reducing redundancy via Theorem 7 of Ref. 2, 
the number of orderings it is necessary to investigate reduces to 360 
for all but 2 of the 22 filters. 

The results of the investigations for all 27 filters will eventually be 
presented. Meanwhile, attention is focused on a typical 13-point filter. 
As an example, a filter with 4 ripples in the passband, 3 ripples in the 
stopband, and passband and stopband tolerances of 0.1 and 0.01, 
respectively, is used. By passband and stopband tolerances is meant 
the maximum height of ripples in the respective frequency bands. 
Figure 5 shows the positions of the zeros of the filter in the upper half 
of the z-plane. Each section of the filter is given a number for identifica- 
tion. The zeros that a section synthesizes are given the same number, 
and these are shown in Fig. 5. Table I shows a list, in order of increas- 
ing noise magnitude, of all 360 orderings investigated and their cor- 
responding output noise variances in units of Q?, computed according 
to Fig. 4. A histogram plot of the noise distribution is shown in Fig. 
6a, and a cumulative distribution plot is shown in Fig. 6b. 

Two characteristics of the histogram shown in Fig. 6a are of special 
importance because they are common to similar plots for all the filters 
investigated. First of all, most significant is the shape of the distribu- 
tion. It is seen that most orderings have very low noise compared to 
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TABLE I—NOISE VARIANCE OF ALL 360 ORDERINGS 
OF A 18-PoIntT FILTER 

















Order Noise Order Noise Order Noise 
263451 1.0983 416253 1.5957 621453 2.4335 
145263 1.1104 341625 1.6081 245613 2.4699 
145362 1.1131 436152 1.6170 134652 2.4854 
163452 1.1382 142635 1.6286 236415 2.5285 
245163 1.1601 412653 1.6354 613425 2.5443 
245361 1.1605 421653 1.6418 314652 2.5643 
362451 1.1834 241635 1.6458 623415 2.5703 
246351 1.2305 243615 1.6524 631425 2.6047 
162453 1.2456 243561 1.6539 632415 2.6073 
361452 1.2561 164352 1.6583 425631 2.6090 
261453 1.2783 264153 1.6817 612435 2.6228 
143652 1.2841 346215 1.6829 621435 2.6461 
146352 1.38245 346125 1.6904 425613 2.6785 
415263 1.3298 364251 1.7040 145632 2.6801 
415862 1.3325 164253 1.7101 134562 2.6991 
243651 1.3356 413562 1.7171 145623 2.6993 
345261 1.3546 342615 1.7177 624351 2.7021 
345162 1.3568 342561 1.7192 326415 2.7167 
246153 1.3652 413625 1.7449 184625 2.7269 
346251 1.3660 431562 1.7483 126453 2.7639 
341652 1.3666 426315 1.7560 314562 2.7779 
425163 1.3687 431625 1.7762 234651 2.7927 
425361 1.3692 412563 1.7835 314625 2.8058 
146253 1.3763 416325 1.7854 634251 2.8110 
163425 1.3797 426135 1.7865 614352 2.8229 
342651 1.4009 364152 1.7869 624153 2.8368 
263415 1.4151 421563 1.7900 614253 2.8747 
142653 1.4160 416235 1.8083 634152 2.8939 
241653 1.4332 412635 1.8480 415632 2.8994 
426351 1.4392 436215 1.8509 216453 2.9085 
346152 1.4489 421635 1.8544 415623 2.9187 
162435 1.4582 436125 1.8585 126435 2.9765 
261435 1.4909 423615 1.8610 324651 2.9809 
361425 1.4976 423561 1.8626 624315 3.0190 
143562 1.4977 264315 1.8638 624135 3.0495 
362415 1.5002 432615 1.8858 614325 3.0644 
413652 1.5034 432561 1.8873 614235 3.0873 
435261 1.5227 264135 1.8943 234615 3.1095 
435162 1.5249 164325 1.8998 234561 3.1110 
143625 1.5256 164235 1.9227 216435 3.1211 
436251 1.5341 364215 2.0208 634215 3.1278 
431652 1.5347 364125 2.0284 634125 3.1354 
416352 1.5439 136452 2.0700 124653 3.2439 
423651 1.5442 316452 2.1489 324615 3.2977 
264351 1.5470 236451 2.2117 324561 3.2992 
246315 1.5474 623451 2.2534 214653 3.3885 
142563 1.5642 632451 2.2904 124563 3.3920 
146325 1.5660 613452 2.3028 124635 3.4565 
432651 1.5690 136425 2.3115 246531 3.4977 
426153 1.5739 631452 2.3632 214563 3.5367 
246135 1.5778 316425 2.3904 246513 3.5672 
341562 1.5803 326451 2.3999 214635 3.6011 
241563 1.5814 245631 2.4004 426531 3.7063 


146235 1.5889 612453 2.4102 426513 3.7758 
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Order 


264531 
264513 
163245 
146532 
146523 
162345 
261345 
361245 
345621 
416532 
345612 
416523 
263145 
164532 
362145 
164523 
435621 
435612 
451263 
451362 
452163 
452361 
453261 
453162 
136245 
624531 
145236 
245136 
316245 
624513 
613245 
415236 
612345 
425136 
631245 
621345 
236145 
145326 
142536 
623145 
241536 
632145 
614532 
614523 
126345 
326145 
415326 
412536 
421536 
345126 
216345 
143526 
245316 
435126 
452631 


Noise 


3.8141 
3.8836 
4.0265 
4.0376 
4.0569 
4.0697 
4.1025 
4.1445 
4.2234 
4.2570 
4.2737 
4.2763 
4.2914 
4.3715 
4.3766 
4.3907 
4.3915 
4.4417 
4.5778 
4.5806 
4.6348 
4.6353 
4.7361 
4.7383 
4.9584 
4.9693 
4.9857 
5.0354 
5.0372 
5.0388 
5.1911 
5.2051 
5.2343 
5.2440 
5.2515 
5.2577 
5.4048 
5.4322 
5.4395 
5.4466 
5.4567 
5.4836 
5.5361 
5.9503 
5.5880 
5.5930 
5.6515 
5.6589 
5.6653 
5.6759 
5.7326 
5.8168 
5.8434 
5.8440 
5.8751 





Order 


341526 
452613 
413526 
345216 
346521 
425316 
431526 
346512 
451632 
451623 
435216 
436521 
436512 
243516 
364521 
342516 
364512 
423516 
432516 
134526 
314526 
124536 
214536 
634521 
634512 
453621 
453612 
234516 
324516 
451236 
452136 
451326 
453126 
452316 
453216 
462351 
463251 
461352 
462153 
463152 
461253 
462315 
462135 
463215 
461325 
463125 
461235 
462531 
462513 
461532 
461523 
143265 
341265 
142365 
241365 


Noise 


5.8993 
5.9446 
6.0362 
6.0375 
6.0426 
6.0520 
6.0674 
6.0928 
6.1475 
6.1668 
6.2055 
6.2106 
6.2609 
6.3368 
6.3805 
6.4021 
6.4307 
6.5454 
6.5701 
7.0181 
7.0970 
7.2674 
7.4120 
7.4875 
7.5377 
7.6049 
7.6551 
7.7939 
7.9820 
8.4532 
8.5101 
8.8996 
9.0573 
9.3181 
9.4189 
11.8333 
11.8896 
11.9658 
11.9680 
11.9725 
12.0176 
12.1501 
12.1806 
12.2064 
12.2073 
12.2140 
12.2302 
14.1004 
14.1699 
14.6789 
14.6982 
16.3035 
16.3860 
16.4329 
16.4501 





Order 


413265 
431265 
463521 
463512 
412365 
421365 
134265 
314265 
243165 
342165 
423165 
432165 
124365 
214365 
234165 
324165 
642351 
643251 
641352 
642153 
643152 
641253 
642315 
642135 
643215 
641325 
643125 
641235 
143256 
341256 
142356 
241356 
413256 
431256 
412356 
421356 
642531 
642513 
134256 
314256 
243156 
342156 
641532 
641523 
423156 
432156 
124356 
214356 
234156 
324156 
643521 
643512 
132645 
312645 
231645 





Noise 


16.5228 
16.5541 
16.5661 
16.6163 
16.6522 
16.6587 
17.5048 
17.5837 
17.7595 
17.8249 
17.9682 
17.9929 
18.2607 
18.4053 
19.2166 
19.4048 
21.7670 
21.8232 
21.8995 
21.9017 
21.9061 
21.9513 
22.0838 
22.1143 
22.1401 
22.1410 
22.1476 
22.1639 
23.0109 
23.0934 
23.1403 
23.1575 
23.2302 
23.2615 
23.3596 
23.3661 
24.0341 
24.1036 
24.2122 
24.2911 
24.4670 
24.5323 
24.6126 
24.6319 
24.6756 
24.7003 
24.9681 
25.1128 
25.9240 
26.1122 
26.4998 
26.5500 
81.4953 
81.5742 
85.2733 
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TasBLE I—Continued. 





Order Noise Order Noise Order Noise 
321645 85.4615 231465 116.6240 465213 180.9850 
123645 87.0142 321465 116.8120 465132 181.3550 
213645 87.1589 123465 118.3650 465123 181.3750 
456231 99.7815 213465 118.5090 465321 182.8770 
456213 99.8510 132456 119.5530 465312 182.9270 
456132 100.2220 312456 119.6320 645231 190.8490 
456123 100.2410 231456 123.3310 645213 190.9180 
456321 101.7430 321456 123.5190 645132 191.2890 
456312 101.7940 123456 125.0720 645123 191.3080 
182465 112.8460 213456 125.2170 645321 192.8110 
312465 112.9250 465231 180.9150 645312 192.8610 





the maximum value possible. In fact, the lowest range of noise vari- 
ance, in this case between zero and 2 Q?, is the most probable range in 
terms of the number of orderings which produce noise variances in this 
range. The distribution is seen to be highly skewed, with an expected 
value very close to the low-noise end, in this case equal to 19.5 Q?. 
In fact, from the cumulative distribution it is seen that approximately 
two-thirds of the orderings have noise variances less than 4 percent 
of the maximum, while nine-tenths of them have noise variances less 
than 14 percent of the maximum. 

The second characteristic is that large gaps occur in the distribution 
so that noise values within the gaps are not produced by any orderings. 
While Fig. 6a shows this effect only for the higher noise values, a more 
detailed plot of the distribution in the range from zero to 28 Q?, as in 
lig. 6c, shows that gaps also occur for lower noise values. Thus noise 
values tend to occur in several levels of clusters. These observations 
provide the general picture of clusters of noise values, the separation of 
which increases rapidly as a function of the magnitude of the noise 
values, thus forming a highly skewed noise distribution. 

The significance of these results is far-reaching. Given a specific 
filter, because of the abundance of orderings which yield almost the 
lowest noise variance possible, it is concluded that it should not be too 
difficult to devise a feasible algorithm which will yield an ordering, 
the noise variance of which is very close to the minimum. Thus, as 
far as designing practical cascade filters is concerned, it really is not 
crucial that the optimum ordering be found. In fact, it may be far 
more advantageous to use a suboptimal method, which can rapidly 
choose an ordering that is satisfactory, than to try to find the optimum. 
The reduction in roundoff noise gained by finding the optimum solu- 
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_Fig. 6—(a) Noise distribution histogram of filter of Fig. 5. (b) Cumulative noise 
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ig. 5. 
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tion is probably, at best, not worth the extra effort from the design 
standpoint. At least up to the present, no efficient method for finding 
an optimum ordering has been found. 

In Section VI, a suboptimal method is presented which, given a filter, 
yields a low-noise ordering efficiently and has been successfully applied 
to a wide range of filters. Before presenting the algorithm, the be- 
havior of roundoff noise with respect to scaling and other filter parame- 
ters will be further investigated. Also, the nature of high-noise and low- 
noise orderings will be discussed, so that they can be more easily 
recognized. 
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Fig. 7—(a) Noise distribution histogram of typical 11-point filter. (b) Noise 
distribution histogram of another 13-point filter. (c) Noise distribution histogram 
of a 15-point filter. 
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Before ending this section, the noise distribution histograms of an 
11-point and one more 13-point filter are presented in Figs. 7a and b. 
These are seen to exhibit all the characteristics discussed above. The 
major difference between the noise distributions of the two 13-point 
filter examples presented lies in the magnitude of the maximum and 
average noise variances. This difference will be accounted for presently. 

Also presented is the noise distribution for a 15-point filter, in Tig. 
7c. The calculation of this distribution involves 2520 different order- 
ings. This histogram shows even stronger emphasis on the distribution 
characteristics discussed and, together with Figs. 6a and 7a, suggests 
that the skewed shape and large-gap properties of the noise distribu- 
tion of a filter become increasingly pronounced as the order of the filter 
increases. Thus it is expected that the results presented can be general- 
ized for higher-order filters. 


3.3 Dependence of Distributions on Transfer Function Parameters 


l’rom all the calculated noise distributions of the previous section, 
the interesting fact is observed that, though different filters may pro- 
duce very different ranges of output noise variances when ordered in 
all possible ways, the noise variances for each filter always distribute 
themselves in essentially the same general pattern. The differences in 
noise variance ranges among different filters is accounted for by investi- 
gating the dependence of noise distributions on parameters which 
specify the transfer function of a filter. 

The noise distributions of several low-pass extraripple filters with 
various values of the parameters, NV, F,, Di, and D2 were computed us- 
ing the methods described. Since all these distributions have the same 
general shape, they can be compared by simply examining their maxi- 
mum, average, and minimum values. A list of all the filters, the noise 
distributions of which have been computed, including those already dis- 
cussed, is presented in Table II. These filters are specified by five pa- 
rameters, namely the four already mentioned, plus N,, the number of 
ripples in the passband. Since all the filters are extraripple filters, it is 
more natural to specify N, than F;. Of course, N, and F are not inde- 
pendent. The maximum, average, and minimum values of the noise dis- 
tributions of each of these filters are listed in Table II. The last column 
in this table will be discussed in Section VI. 

Filters numbered 1 to 5 in Table II are very similar except for their 
order in that they all have identical passband and stopband tolerances 
and approximately the same low-pass bandwidth. The maximum, 
average, and minimum values of their noise distributions are plotted 
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TABLE [I—LisT oF FILTERS AND THEIR NOISE 
DISTRIBUTION STATISTICS 


















































Noise Variance 

# N N> Py dD, D, Max | Avg Min | Algl 
1 7 2 0.212 | 0.1 0.01 1.24 0.84 0.37 0.37 
2 9 3 0.281 | 0.1 0.01 6.26 2.54 0.73 0.73 
3 11 3 0.235 | 0.1 0.01 19.41 4.79 0.68 0.68 
4 13 4 0.279 | 0.1 0.01 192.86 | 19.55 1.10 1.10 
5 15 4 0.244 | 0.1 0.01 923.63 | 54.45 1.02 1.16 
6 13 3 0.100 | 0.001 | 0.001 15.84 3.01 0.65 0.69 
7 13 4 0.261 | 0.05 0.004 | 119.48} 12.91 0.96 1.02 
8 13 1 0.012 | 0.01 0.01 9.91 1.61 0.32 0.35 
9 13 2 0.067 0.01 0.01 16.30 2.94 0.44 0.47 

10 13 3 0.138 | 0.01 0.01 42.63 5.94 0.71 0.73 

11 13 4 0.213 | 0.01 0.01 69.76 8.52 0.82 0.91 

12 13 5 0.288 | 0.01 0.01 76.43 | 11.01 1.44 1.52 

13 13 6 0.364 | 0.01 0.01 52.54 | 10.33 1.92 2.43 

14 13 3 0.201 | 0.1 0.01 96.25 | 12.09 0.81 

15 13 3 0.179 | 0.05 0.01 69.26 9.02 0.76 

16 13 3 0.154 0.02 0.01 50.63 6.87 0.72 

17 13 3 0.123 | 0.005 | 0.01 37.36 5.33 0.70 

18 13 3 0.106 | 0.002 | 0.01 32.83 4,80 0.69 

19 13 3 0.095 | 0.001 | 0.01 30.53 4.53 0.69 

20 13 3 0.124 | 0.01 0.1 132.57 | 17.56 1.02 

21 13 3 0.129 | 0.01 0.05 85.84 | 11.45 0.83 

22 13 3 0.135 | 0.01 0.02 54.94 7.47 0.75 

23 13 3 0.141 | 0.01 0.005 35.59 5.07 0.68 

24 13 3 0.144 | 0.01 0.002 26.44 4.37 0.68 

25 13 3 0.146 | 0.01 0.001 22.52 4.07 0.70 

26 15 4 0.185 | 0.01 0.01 417.08 | 27.38 1.00 

27 15 4 0.255 | 0.1 0.001 | 601.83 | 35.15 1.02 














on semilog coordinates in I*ig. 8a. It is seen that all these statistics of 
the distributions have an essentially exponential dependence on filter 
length. The less regular behavior of the minimum values is believed 
to be caused by differences in bandwidth (/1) among the filters. 
Tigure 8b shows a similar plot of the same distribution statistics for 
filters numbered 8 to 13 as a function of F;. These filters have identical 
values of N, Di, and De, and represent all six possible extraripple 
filters that have these parameter specifications. From Fig. 8b it is 
seen that with those parameters mentioned above held fixed, the noise 
output of a cascade filter tends to increase with increasing bandwidth. 
Filters numbered 14 to 25 all have fixed values of N, N,, and either 
D, or Dz. Plots of the distribution statistics of these filters as func- 
tions of D, and D, are shown respectively in Figs. 8c and 8d. These 
plots indicate that, as the transfer function approximation error for a 
filter decreases, so does its noise output. Though the plots are made 


ROUNDOFF NOISE IN FIR FILTERS 363 


holding NV, rather than F’ fixed, it is seen that, at least for the filters 
used in Fig. 8d, bandwidth increases with decreasing approximation 
error. Since the noise output of a filter is found to increase with band- 
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Fig. 8—(a) Output noise variance as a function of filter length. (b) Output noise 
variance as a function of bandwidth. (c) Output noise variance as a function of pass- 
band approximation error. (d) Output noise variance as a function of stopband 
approximation error. 
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width, it is expected that noise would still decrease with stopband 
tolerance D» if F, were fixed instead of N,. In any event, the variation 
of F, among these filters is small. 

Figures 8b to 8d are all plots of statistics for 13-point filters. Notice 
how the maximum, average, and minimum curves all tend to move 
together. In particular, the average curve almost always stays approxi- 
mately halfway on the logarithmic scale between the maximum and 
minimum curves. This phenomenon is, of course, simply a manifesta- 
tion of the empirical finding that noise distributions of different filters 
have essentially the same shape independent of differences in transfer 
characteristics. 

To summarize, it has been found experimentally that with other 
parameters fixed, the roundoff noise output of a filter tends to increase 
with increases in all four independent parameters N, F';, Di, and De 
which specify its transfer function. In particular, noise output tends to 
grow exponentially with N. It was not shown that the noise output of 
a filter with a fixed ordering and scaled a given way always varies in 
the way indicated when its transfer function parameters are perturbed. 
What has been shown is perhaps a more useful result from the design 
viewpoint. These findings imply that, other things being equal, a 
transfer function with, for instance, a higher value of Dz, is likely, 
when realized in a cascade form, to result in a higher noise output than 
a transfer function with a smaller value of D» realized by the same 
method. Though these results were found using only low-order filters, 
it is expected they could be generalized for higher-order filters as well. 
Section VI will present experimental evidence to confirm this 
expectation. 


IV. COMPARISON OF SUM SCALING AND PEAK SCALING 


The claim was made earlier that the results obtained on the noise 
distribution of filters ought to be quite independent of whether sum 
scaling or peak scaling is used. This claim will now be sustained 
heuristically and experimentally. 

Let a denote the ratio of the maximum gain (over all frequencies) 
of a low-pass extraripple filter scaled by sum scaling to the maximum 
gain of the same filter peak scaled. Then it must be true that a S 1, 
since by definition the maximum gain for peak scaling is exactly one, 
while for sum scaling it must be no more than one if class 2 signals 
(which are a subset of class 1) are to be properly constrained (by 
Theorem 3, Ref. 2). Furthermore, it can be easily shown! that 
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a = 1 — 2e where e« is defined as 
>» r(k) 
k 


¢ = ————_ 3 
dX I A(k)| o 


with {h(k)} being the impulse response of the filter and {r(k)} being 
the magnitude of the negative portion of {h(k)}, 1e., 


—h(k) h(k) < 0 


anime, h(k) = 0 


(4) 


For a low-pass filter, the envelope of {h(k)} has the general shape 
of a truncated (sin x) /x curve, hence ¢ is expected to be a small number. 


TABLE [][—List oF FILTERS AND THE RESULTS 
OF ORDERING ALGORITHMS 








Noise Variance 


























D1 = 0.01 . Sum 

D2 = 0.001 Peak Scaling Scaling 
# N N> Fy a Alg 1 Alg 2 Alg 1 
23 | 13 4 | 0.219 0.65 1.25 1.26 0.90 
29 | 15 4 | 0.193 0.68 123 1.22 1.02 
30 | 17 5 | 0.230 0.61 1.99 2/49 137 
31 | 19 5 | 0.207 0.64 1.93 1.92 147 
32 | 21 6 | 0.236 0.59 250 2.61 1.58 
33 | 23 6 | 0.216 0.61 2.57 2.91 1.77 
34 | 25 7 | 0.240 0.57 3.75 3.62 2°35 
35 | 27 7 | 0.223 0.59 3.95 4.11 2.45 
36 | 29 8 | 0.243 0.55 4.54 5.04 2.67 
37 | 31 g | 0.297 0.57 5.27 5.88 274 
38 | 33 9 | 0.244 0.54 781 6.67 4°59 
39 | 35 9 | 0.231 0.55 6.01 6.43 3.72 
40 | 33 1 | 0.005 1.0 0.47 0.48 0.53 
41 | 33 2 | 0.029 0.82 0.60 0.67 0.60 
42 | 33 3 | 0.059 0.73 0.89 1.00 0.80 
43 | 33 4 | 0.090 0.68 1.43 1.36 1.16 
44 | 33 5 | 0121 0.63 2.29 1.84 71 
45 | 33 6 | 0.152 0.60 2.48 2:70 1.61 
46 | 33 7 | 0.183 0.58 3.47 3.37 230 
47 | 33 gs | 0.214 0.61 4.72 5.23 3.38 
4g | 33 | 10 | 0.275 0.52 10.04 8.16 4.83 
49 | 33 | 1 | 0.305 0.52 15.68 11:35 8:30 
50 | 33 | 12 | 0.334 0.50 13.43 14.88 6.27 
51 | 33 13 | 0.363 0.50 21.35 17.62 9.14 
52 | 33 | 14 | 0.392 0.50 41.64 | 31.41 15.40 
53 | 33 | 15 | 0.419 0.51 55.20 | 41.13 22.12 
54 | 33 | 16 | 0.448 0.53 89.52 65.66 38.23 
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TABLE [V—List or FILTERS AND THE RESULTS 
OF ORDERING ALGORITHMS 








Noise Variance 


N = 33 . Sum 
N,=8 Peak Scaling Scaling 
# Fy dD, D, a Alg 1 Alg 2 Alg 1 
55 0.211 0.01 0.002 0.59 5.63 5.33 3.34 
56 0.208 0.01 0.005 0.57 5.13 5.69 3.18 
57 0.205 0.01 0.01 0.56 5.05 5.27 3.34 
58 0.202 0.01 0.02 0.55 7.63 8.31 4.01 
59 0.197 0.01 0.05 0.53 11.34 12.53 6.92 
60 0.193 0.01 0.1 0.51 46.33 22.99 16.88 
61 0.238 0.1 0.01 0.58 9.90 9.01 5.61 
62 0.227 0.05 0.01 0.58 8.91 7.35 5.52 
63 0.214 0.02 0.01 0.56 8.87 5.75 4.32 
64 0.196 0.005 0.01 0.56 5.47 4.69 3.68 
65 0.185 0.002 0.01 0.56 5.95 4.08 3.41 
66 0.178 0.001 0.01 0.57 4.10 4.11 2.85 


In fact, it is easily shown! that, for a low-pass filter, if the passband 
tolerance D; is much less than the maximum passband gain, as is 
usually the case, then to an excellent approximation a = 1 — 2e. It is 
now shown that if o? is the output noise variance of a filter with sum 
scaling and o” is the output noise variance of the same filter except 
with peak scaling, then o? = a?(c”). To show this, the optimality 
properties for sum scaling and peak scaling, proved in Ref. 2, Theorem 
4, are invoked. Since the gain of the sum-scaled filter is a times that of 
the peak-scaled filter, the ratio of their signal-to-noise ratios for class 
2 inputs must equal a times the inverse ratio of their rms noise values 
(square root of variance), i.e., with S/N for sum scaling and S/N’ 
for peak scaling, 
S/N a’ 


S/N’! ©) 





But since peak scaling is optimum for class 2 inputs and class 2 is a sub- 
set of class 1, then S/N’ = S/N. Thus o? 2 a’o”. For an alternative 
derivation see Ref. 1. 

In Tables III and IV, a list of filters and some results of Section VI 
are presented. Together with these results are listed measured values 
of a for each filter. Observe that, for these typical filters, a ranges from 
0.5 to 1. Furthermore, for each filter, the last and third to last columns 
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N =13 


F, = 0.201 


F, = 0.301 
D, =0.1 


D, = 0.01 


NOISE VARIANCE FOR PEAK SCALING 





0 10 20 30 40 50 60 70 80 90 100 
NOISE VARIANCE FOR SUM SCALING 


Fig. 9—Peak scaling versus sum scaling noise output comparison for a typical filter. 


of Tables III and IV list the noise variances that result from the same 
ordering using sum scaling and peak scaling respectively. Comparing 
these, it is seen that, in almost every case, o? S (o”).? In particular, 
this is true if a is not too close to 1.0. The case where o? > (c”) and 
a = 1 is filter number 40 in Table III. However, except for the un- 
interesting cases of filters with all zeros on the unit circle, in general, 
a < 1, and it is expected that o? S$ (o”). Thus for all practical pur- 
poses it can be assumed that 

2 


q 


a 


IA 
| 
IA 


— $1. (6) 


to 


From eq. (6), it is seen that the output noise variance for a filter 
with sum scaling is comparable, at least in order of magnitude, to that 
for the same filter ordered the same way with peak scaling applied. 
In fact, experimental results show that given a filter, the noise vari- 
ances for sum scaling and peak scaling are in an approximately con- 
stant ratio for almost all orderings. An example of this result is shown 
in Fig. 9, where the noise variances for sum scaling and peak scaling of a 
typical filter are plotted against each other for each ordering. The 
resulting points are seen to form almost a straight line with slope 
approximately equal to 2, so that essentially o’? = 20? for all orderings 
of this filter. 

Thus the noise distributions of the previous section are essentially 
unchanged if peak scaling is used instead of sum scaling. To illustrate 
this, Fig. 10 shows the noise distribution for the filter of Fig. 6 with 
peak scaling used instead of sum scaling. 
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TOTAL ORDERINGS = 360 
W = 5.05 


NUMBER OF ORDERINGS 
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Fig. 10—Noise distribution histogram of filter of Fig. 5 using peak scaling. 


The evaluation of noise variances with peak scaling is done in essen- 
tially the same way as that described in Fig. 4. Using a 128-point FFT 
to evaluate two transforms at a time (exploiting real and imaginary 
part symmetries) to give the maxima of the F;(e*) for 360 orderings, 
the computations for peak scaling were found to require four times as 
much time as that for sum scaling. 


Vv. AN INTUITIVE EXPLANATION OF ROUNDOFF NOISE DEPENDENCE ON 
ORDERING IN TERMS OF SPECTRAL PEAKING 


That roundoff noise is distributed in the way shown with respect to 
orderings for a filter is an intriguing fact which is by no means obvious. 
The dependence of roundoff noise on ordering involves complicated 
matters like differing spectral shapes of different combinations of 
individual filter sections and the interactive scaling of signal levels 
within a filter necessitated by dynamic range limitations. As such, this 
dependence is much too complicated to visualize intuitively. It is 
proposed that the relative level of roundoff noise in a filter is adequately 
determined in order of magnitude by the amount of peaking in certain 
subfilter spectra. Thus it will be shown that, since the dependence on 
ordering of the amount of peaking of these spectra is not too difficult 
to visualize, by judging the relative amount of peaking of these spectra, 
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the relative merit of an ordering in terms of high-noise or low-noise 
output can be determined by inspection. These findings explain the 
general shape of noise distributions. 

Given a linear phase filter with z-transform H(z), define transfer 
functions H;(z),7 = 1, ---, Ns, to have the property that H;(z) is pro- 
portional to the transfer function for the 7th section of the filter, and 
each | H;(e%)| for all i has a maximum over w equal to C, where C is 
chosen so that the overall filter frequency response H (e*”) has a maxi- 
mum in magnitude equal to one. Clearly C = 1. Define transfer 
functions 


Ade) = TL Hi) 

e (7) 
Bie) = I Hy) 

jeitl 


and define a number Pk to be the largest value of max.,|.A;(e%)| or 
max, | B;(e%)| for all 7. It will be argued that, given an ordering, a 
large value of Pk indicates a high-noise output, while a low value of 
Pk indicates a low-noise output. 

To see this, define G,(e*) to be B;(e”) with its maximum in mag- 
nitude over w normalized to unity. Then it can be easily shown! that if 


A;= max | A,(e%) | 


B; = max | B.(e%) | 


Q? 1 27 = 
= | | G;(e%) | *dw 
12 Qn 0 


(8) 


(where k; is the number of noise sources in the 7th section), then the 
output noise variance due to the 7th section is given by 


ABO, ee pen esis (9) 


For the moment assume that C; is a constant factor independent of 
ordering. Then o7 is proportional to (A;B,)?. Note that for any 7, A; 
and B; are the maxima of two functions the product of which is H (e%). 
Furthermore, for some 7 either A; = Pk or B; = Pk. Now suppose 
Pk > C. Without loss of generality, it may be assumed A; = Pk. 
Then argue that A;B; > C. 

Clearly A; = |A;(e%o)| for some wo. Now A;(z)B;:(z) = H(z), and 
H(z) is a function with zeros only in the z-plane other than the origin. 
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Also, at least in the case of well-designed band-select filters, the zeros 
of H(z) are well spaced and spread out around the unit circle. The 
zeros of a typical filter are shown in Fig. lla. Furthermore, |H (e%*) | 
<1. Thus in order for |A;(e)| to have a large peak at wo, several 
zeros of H(z) which occur in the vicinity of 2 = e+’. must be missing 
from A;(z), while most of the remaining zeros must be in A;(z). This 
means that B;(z) has a concentration of zeros around e+4#o. By the re- 
sult of Theorem 6(z) in Ref. 2, which says that the maximum of the 
magnitude frequency response of a filter section occurs at either w = 0 
or w = x depending on whether the zeros synthesized are in the left half 
or the right half of the z-plane, it is seen that most factors of B;(e%) 
must have maxima in magnitude which occur at exactly the same w. 


IMAGINARY PART OF 2 





IMAGINARY PART OF Z 


—2.0 —1.6 —1.2 —0.8 —0.4 0 0.4 0.8 1.2 1.6 2.0 
REAL PART OF Z 


Fig. 11—(a) Zero positions of a typical 67-point filter. (b) Zero positions of filter 
number 14 of Table IT. 
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Now C is found to be an increasing function of Ni, where typically 
C > 2. Hence |B;(e*) | is very likely to have a peak which is at least 
1, or B; 2 1. Thus A;B; > C. By the same token, if B; = Pk and 
Pk > C, then A;:B; > C. 

Hence if Pk >> C, then for at least one 7, o% = (A,B;)?C; where 
A:B;>>C. Compared with a nominal value of say A;B; = C, the 
resulting difference in output noise variance can be great. When Pk 
takes on its lowest possible value, viz., Pk = C, the o%’s are compara- 
tively small for all 2, hence it is expected that the resulting o? is among 
the lowest values possible. Thus there exists a correlation between 
high values of Pk and high noise, and low values of Pk and low noise. 

Concerning the assumption that C; is constant independent of order- 
ing, it is reasonable as long as only order-of-magnitude estimates are 
of interest. Since by definition max.|G;(e%)| = 1 independent of 
ordering and 7, it can be expected that variations in C; with ordering 
are much less than variations in (A;B;). 

Based on these results, it can be concluded that an ordering which 
groups together, either at the beginning or at the end of a filter, several 
zeros all from either the left half or the right half of the z-plane is 
likely to yield very high noise. This observation is based on the fact 
that since zeros from the same half of the z-plane produce frequency 
spectra the maxima of which occur at exactly the same w, several zeros 
from the same half of the z-plane can build up a large peak in the pro- 
duct of their spectra H;(e*). On the other hand, a scheme which 
orders sections so that the angle of the zeros synthesized by each 
section lies closest to the w at which the maximum of the spectrum 
of the combination of the preceding sections occurs is likely to yield 
a low-noise filter. 

The above observations are found to be true for all the filters the 
noise distributions of which were investigated. For example, from 
Table I, it is seen that those orderings which group together all three 
sections 4, 5, and 6 of the filter of Fig. 5 either at the beginning or at the 
end of the filter are precisely those which have the highest noise. 
Namely, they account for all noise variance values above 26.6 Q?. 
Using the results on the noise distribution of a filter and the results 
of this section, it can thus be said that the comparatively few orderings 
of a filter which have unusually high noise can be avoided simply by 
judiciously choosing zeros for each section so that no large peaking 
in the spectrum, either as seen from the input to each section or from 
each section to the output, is allowed to occur. In particular, this can 
be done by ensuring that from the input to each section the zeros 
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synthesized well represent all values of w, i.e., the variation in the 
density over w of zeros chosen should be minimal. 

These observations have the implication that low-noise orderings 
for a filter are those which choose zeros in such an order that they 
‘Jump around”’ about the unit circle and are well “‘interlaced,’’ whereas 
high-noise orderings are those which allow clusters of adjacent zeros 
to be sequenced adjacently in the filter cascade. Since there are cer- 
tainly far more ways to sequence zeros so that they satisfy the former 
property than the latter, it is clear why most orderings of a filter have 
low noise. 

The values of Pk for all orderings of several filters were measured,! 
and they show good correlation with o?, thus supporting these argu- 
ments. For reasons of space, the results are not tabulated here. How- 
ever, Figs. 12 and 13 show plots of the spectra from each section to the 
output of a high-noise and a low-noise ordering for a typical 13-point 
filter, the zeros of which are shown in Fig. 11b. Both orderings are 
peak scaled, so that the spectra from the input to each section of the 
filter have maxima equal to one. Thus, in reference to eq. (9), A:B; is 
equal to the maximum of the spectrum from section (7 + 1) to the 
output. The ordering of Fig. 12 has o? = 186 Q?, while that of Fig. 13 
has o? = 1.1 Q*. It is seen that, as expected, A,B; has a large value for 
at least one 7 in the high-noise ordering, reaching a value of 60, while 
for the low-noise ordering A;B; < 2.2 for all 7. Furthermore C;, which 
is proportional to the integral of the square of the spectrum from section 
(¢ + 1) to output with its maximum normalized to unity, does not 
vary too much between the two orderings. Finally, note that the spec- 
trum of each section in the low-noise ordering does indeed tend to 
suppress the peak in the spectrum of the combination of previous 
sections. Thus the arguments of this section are supported. 


VI. AN ALGORITHM FOR OBTAINING A LOW-NOISE ORDERING FOR THE 
CASCADE FORM 


An extensive analysis of roundoff noise in cascade form FIR filters 
has been presented in this paper and in Ref. 2. However, an investi- 
gation of roundoff noise would not be complete without studying the 
practical question which in the first place had motivated all the 
analyses and experimentation. The question is, given an FIR transfer 
function desired to be realized in cascade form, how does one sys- 
tematically choose an ordering for the filter sections so that roundoff 
noise can be kept to a minimum? 

A partial answer to this question has already been given in the 
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FILTER SPECTRUM FROM SECTION 4 TO OUTPUT 
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FILTER SPECTRUM FROM SECTION 5 TO OUTPUT 
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FILTER SPECTRUM FROM SECTION 6 TO OUTPUT 
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Fig. 12—Series of spectra from section 7 to the output where z = 2, 3, 4, 5, 6, for a 
high-noise ordering. 
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Fig. 13—Series of spectra from section 7 to the output where z = 2, 3, 4, 5, 6, for a 
low-noise ordering. 
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previous section. However, no completely systematic method has yet 
been devised for selecting an ordering for a filter guaranteed to have 
low noise. Ultimately, one wishes to find an algorithm which, when 
implemented on a computer, can automatically choose a proper order- 
ing in a feasible length of time. 

Avenhaus has studied an analogous problem for cascade IIR filters 
and has presented an algorithm for finding a “favorable” ordering of 
filter sections. His algorithm consists of two major steps; a ‘‘pre- 
liminary determination” and a ‘‘final determination.”’ In this section 
an algorithm is described for ordering FIR filters which is based upon 
the procedure used in the “preliminary determination” step of Aven- 
haus’ algorithm. It has been found that a procedure appended to the 
proposed algorithm similar to Avenhaus’ final determination step 
adds little that is really worth the extra computation time to the 
already very good solution obtainable by the first step. Hence such a 
procedure is not included in this algorithm. 

No statement was made by Avenhaus as to what range of noise 
values can be expected of filters ordered by his algorithm, nor did he 
claim that his algorithm always yields a low-noise ordering (relatively 
speaking, of course). However, based on the results of Sections III 
through V, it will be argued heuristically that the proposed algorithm 
always yields filters which have output noise variances among the 
lowest possible. Together with extensive experimental confirmation, 
these arguments provide confidence that the proposed algorithm 
produces solutions that are very close to the optimum. 

Application of Avenhaus’ procedure to FIR filters allows the intro- 
duction of modifications which reduce significantly the amount of 
computation time required. Also, while IIR filters seldom require 
a higher order than the classic 22nd-order bandstop filter quoted by 
Avenhaus, practical FIR filters can easily require orders over 100. 
Though the same basic algorithm should still work for high orders, care 
must be exercised in performing details to avoid large roundoff errors in 
the computations. Through proper initialization, the proposed algo- 
rithm has been successfully tested for filters of order up to 128. 


6.1 Description and Discussion of Basic Algorithm 


The basic procedure or algorithm proposed by Avenhaus is simply 
the following. To order a filter of N, sections, begin with 7 = N, and 
permanently build into position 2 in the cascade the filter section which, 
together with all the sections already built in, results in the smallest 
possible variance for the output noise component due to noise sources 
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in the zth section of the cascade. Because in the FIR cascade form noise 
is injected only into the output of each section, for FIR filters the 
procedure needs to be modified by considering the output noise due to 
the section in position 7 — 1 rather than 7 when choosing a section for 
position 7. But the 7th section is determined before the (7 — 1)th 
section, hence the number of noise sources at the output of the 
(¢ — 1)th section is unknown at the time that a section for position 
tis to be chosen. This problem is overcome by assuming all sections 
to have the same number of noise sources. Then oc? is simply propor- 
tional to >, 9%(k) independent of what the 7th section is, where 
{g:(k)} is the impulse response of the part of the filter from section 
(¢ + 1) to the output. 

Hence the revised basic algorithm for ordering FIR cascade filters is: 
beginning with 1 = N,, permanently build into position 7 the section 
which, together with the sections already built in, causes the smallest 
possible value for >> x. 9?—1(k). Once this basic algorithm is determined, it 
is necessary only to decide on a scaling method and a computational 
algorithm for accomplishing the desired scaling and noise evaluation 
before an ordering algorithm is completed. Prior to discussing these 
issues, let us consider why the basic algorithm described above is al- 
ways able to find a low-noise ordering. 

The reason why the algorithm might not be able to find a low-noise 
ordering is that rather than minimizing > oc? directly, it minimizes 
each oj individually where for oj, 1 S$ 7 SN, —1, the search is 
essentially conducted over only (7+ 1)! out of the total of N,! 
possible orderings. Now this set of (7 + 1)! orderings depends on which 
sections were chosen for positions 7 + 2 to N, in the cascade if 7 < N, 
— 1. Hence in choosing a section for position 7, previous choices might 
prevent attainment of a sufficiently small value for o7_. 

The basis for the following arguments is presented in Section V. 
Let H(z) be an appropriately scaled filter. Given 1,1 S15 N,—1, 
suppose o; is small for all 7 = J. Then the zeros of [[*,4; H(z) must 
be well spread around the unit circle in the z-plane since a clustering 
would cause large peaking in [[*,4, H.(e) for some k = 1, hence 
a large value of o7. But this means that the remaining zeros of H(z), 
namely those in []‘~, H;(z), must also be well spread around the unit 
circle, since the zeros of H(z) are distributed almost uniformly around 
the unit circle. Hence it ought certainly to be possible to find some pair 
of zeros in [[$_, H;(z) which, when assigned to position 1, causes little 
peaking in [[j=} A; (e%) or [[*, H,(e%), and thus results in a small 
value for oj_,. By induction, then, o% can be chosen small for all 7. 
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For small / it is true that there are very few zeros left as candidates 
for position 1, but in these positions little peaking in the spectra can 
occur since the overall spectrum []*; H;(e”) must be a well-behaved 
filter characteristic. Typically in a high-noise ordering, o7 reaches a 
peak for 1 somewhere in the middle between 1 and N,, while oj for 
small J has little contribution to o?, the total output noise variance. 
Hence the choice of sections for small J is not too crucial. Of course, the 
eligible candidates are still well-spaced zeros as for larger J, so that 
peaking should not be a problem. 

Note that the reason the algorithm works so well is tied in with the 
result of Section III that most orderings of a filter have comparatively 
low noise. Because it is not difficult to find low-noise arrangements of 
zeros, >. a7 can be minimized approximately by minimizing each o% 
independently, searching over a much smaller domain. If the sum 
> o% could not be segmented, searching for a minimum would be 
essentially an impossible task because of time limitations. 


6.2 Two Versions of the Complete Algorithm 


Having discussed why the basic algorithm works, the practical 
problem of implementing it is now discussed. First of all, there is the 
choice of scaling method to use in computing the +; g7(k). As in the 
calculation of noise distributions in Section III, sum scaling is to be 
preferred since it can be carried out the fastest. Figure 14 shows a 
flow chart of the ordering algorithm in which sum scaling is employed. 
Calculation of o? (Nx in the flow chart) 1s done exactly the same way 
as in the algorithm of Fig. 4. 

Using this ordering algorithm, over 50 filters have been ordered and 
the noise variances in units of Q? (Q = quantization step size) of the 
resulting filters are shown in the last columns of Tables II through 
IV. Note that these noise variances are computed with sum scaling 
applied to the filters. The corresponding noise variance values for peak 
scaling have also been computed for the filters of Tables III and IV. 
These are shown in the third to the last column of the tables. The 
comparability of these noise values to those for sum scaling is a con- 
firmation of the results of Section IV. 

For an alternative implementation of the basic ordering algorithm, 
peak scaling can be used. To distinguish between the two different 
resulting algorithms, the former (sum scaling) will be referred to as 
alg. 1 and the latter as alg. 2. The only change needed in Fig. 14 to 
realize alg. 2 rather than alg. 1 is to replace >| f:(k)| by max. | F;(e*) | 
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for given 7 whenever it appears. Results of using alg. 2 on the filters of 
Tables III and IV are shown in the next to last column of those tables. 
Observe that though the two algorithms in general yield different 
orderings for a given filter, the resulting noise variances are comparable. 
Thus, using both alg. 1 and alg. 2, two separate low-noise orderings 
for a given filter can be obtained. Further discussion of the results 
will be given shortly. 

Even with a scaling method decided upon, the questions still remain 
of how +, g2(k) and | fi(k)| or max,|F;(e”)| are to be computed 
and how the sequence {IO(z)} which describes how the filter sections 
are ordered (see Fig. 14) is to be initialized. In obtaining the results of 
Tables III and IV, the following has been done: >; g3(k) and | f:(k) | 
were computed by evaluating {g;(k)} or {f,;(k)} through simulation in 
the time domain (i.e., convolution) and max,,| F;(e%) | was determined 
by transforming {f;(k)} via an FFT and then finding the maximum 
value. Finally, {10 (2)} was initialized to IO(2) = 7,7 = 1, ---, Ns. It 
will be seen that these procedures must be modified for higher-order 
filters. But meanwhile, the implications of these procedures in terms of 
dependence of computation time on filter length is considered. 

Clearly, in algorithmically computing the impulse response of an 
N-point filter via convolution, the number of multiplies and adds 
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Fig. 15—Computation time versus filter length for ordering algorithm. 
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required to calculate each point varies as N, hence the time required 
to evaluate the entire impulse response must vary approximately as 
N?. Now in the basic algorithm there are two nested loops, where the 
number of times the operations within the inner loop are performed 
is given by 

Ns i 

See NAN, + 1) 
ind 2 


N2 


—~ g 





Clearly for alg. 1, the evaluation of >°x|fii1(k)| and Dox g?_1(k) 
dominates all operations within the inner loop in terms of time re- 
quired. Since the total number of points that must be evaluated in 
order to compute {fi-1(k)} and {g:1(k)} together turns out to be a 
constant independent of 1, the combined operations must have approxi- 
mately an N? time dependence. Hence it is predicted that the com- 
putation time required for alg. 1 must be approximately proportional 
to N*. This prediction is verified in Fig. 15, where computation time 
for alg. 1 on the Honeywell 6070 computer is plotted against N on 
log-log coordinates for various values of N. As expected, these points 
lie on a straight line with a slope very nearly equal to 4. 

For alg. 2, exactly the same procedures as in alg. 1 are carried out 
except that after each evaluation of {f;(k)} an FFT is performed. 
Thus for a given N, alg. 2 always requires more time than alg. 1, with 
the exact difference depending on the number of points employed in 
the FFT. 


6.3 Modification of Algorithm for Higher Filter Orders 


For filters of length greater than approximately 41, it is found that 
accuracy in the evaluation of impulse response samples by the methods 
described rapidly breaks down. This phenomenon is chiefly due to the 
fact that the initial ordering used is a very bad one. In particular, it is 
not difficult to see that this ordering (i.e., IO(¢) = 7) has a noise 
variance which is among the highest possible and which increases at 
least exponentially with N. Thus all attempts at evaluating the impulse 
response of the filter by simulation in the time domain are marred by 
roundoff noise. 

A natural possibility for resolving this problem is to perform cal- 
culations in the frequency domain. This has been tried as a modifica- 
tion to alg. 2. In particular, rather than computing F;(e”) from 
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{fi(k)}, it is evaluated as a product of H,(e%), 1 = 1, ---, 7, where 
Hi(e*”) is computed from the coefficients of section / via an FFT. In 
this way the accuracy problem was solved, but computation time in- 
creased significantly. As an example, the 67-point filter listed in Table 
V was ordered using this method. The resulting noise variance was a 
reasonable 26.6 Q?, but even with a 256-point FFT the computation 
time required amounted to 7.2 minutes, more than 7 times that re- 
quired for alg. 1 to order the same filter. 

A far better solution is as follows. The conclusion of Section IJI-that 
most orderings of a filter have relatively low noise—means that, if an 
ordering were chosen at random, it ought to have relatively low noise. 
The strategy is then to use a random ordering as an initial ordering 
for alg. 1. A given ordering of a sequence of numbers {IO(7z),7 = 1, ---, 
N,} can be easily randomized using the following shuffling algorithm :" 


Step 1: Set 7<N,. 

Step 2: Generate a random number U, uniformly distributed be- 
tween zero and one. 

Step 3: Set k[jU]+ 1. (Now k is a random integer between 1 
and j.) Exchange IO(k)  IO()). 

Step 4: Decrease 7 by 1. If 7 > 1, return to Step 2. 


By adding a step to randomize the initial ordering IO(z) = 7 in alg. 
1, the inaccuracy problem was eliminated. The interesting question 
now arises that, since most orderings of a filter have relatively low 
noise, can we not obtain a good ordering simply by choosing one at 
random? The answer is yes in a relative sense, but, as will shortly be 
seen, a random ordering is far from being as good as one which can be 
obtained using the ordering algorithm. 

The extra step of randomizing the initial ordering for alg. 1 requires 
negligible additional computation time, and a filter with impulse re- 
sponse length as high as 129 has been successfully ordered in this way. 
The time required to order this filter was approximately 13.5 minutes. 
Except for time limitations, there is no reason why even higher-order 
filters cannot be similarly ordered. Further results are described in the 
next section. 


6.4 Discussion of Results 


Note in Table II that alg. 1 can result in a noise variance which is 
very close to the minimum, if not the minimum. From this observation 
and the conclusions of Section III on the dependence of the minimum 
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noise variance for a filter on different parameters, we can be quite con- 
fident that the noise variances shown in Tables III and IV are also very 
close to the minimum possible. The filters of Tables III and IV were 
chosen intentionally to show the dependence of the results of the order- 
ing algorithms on various transfer function parameters. It is seen that 
the noise variances indeed behave in the way that would be expected 
from the results of Section III. In particular, c? is seen to be essentially 
an increasing function of N, F;, Di, as well as Ds. The results of 
Tables III and IV are then a confirmation of the expectation that 
the conclusions of Section III on the general dependence of noise on 
transfer function parameters can be generalized to higher-order filters. 

The results of using the modified alg. 1 (denoted alg. 1’) on a few 
filters are shown in Table V. Also shown in this table, for comparison, 
are the noise variances of these filters when they are in the sequential 
ordering [O(z) = 7 (where computable within the numerical range of 
the computer) as well as when they are in a random ordering (obtained 
by randomizing {IO(z)} where IO(¢) = 7, as described above). Be- 
cause of the potentially very large roundoff noise encounterable in 
these orderings, the noise variances were computed using frequency 
domain techniques. In particular, each H;(e”) is evaluated via an 
FFT, peak scaling is then performed, and finally o7 is computed via 
1/(2r) S?™| G;(e%) |2dw rather than +; 93(k). 

From Table V it is seen that though the noise variances of the 
random orderings are certainly a great deal lower than those of the 
corresponding sequential orderings, they are far from being as low as 
those obtained by alg. 1’. Thus it is certainly advantageous to use alg. 1 
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(or 1’) to find proper orderings for filters in cascade form. In all the 
examples given in Tables II through V, one can do little better in 
trying to find orderings with lower noise. With the possible exception 
of the uninteresting wideband filter, number 54 in Table III, all 
filters have less than 4 bits of noise after ordering by alg. 1, while the 
great majority have less than 3 bits. Thus it is not expected. that these 
noise figures can be further reduced by much more than a bit or so. 

Finally, in practice, cascade FIR filters of orders over approximately 
50 are generally of little interest since there exist more efficient ways 
than the cascade form to implement filters of higher orders. For filters 
of, at most, 50th order, the computation time required for alg. 1 is less 
than 20 seconds on the Honeywell 6070 computer. Thus alg. 1 (or 1’) 
is a very efficient means for ordering cascade filters. 


VII. SUMMARY 


In this paper, experimental results have been presented which show 
that most orderings for an FIR filter in cascade form have very low 
noise relatively, and that the shape of the distribution of noise with 
respect. to ordering is essentially independent of transfer function 
parameters as well as method of scaling (sum or peak). An explanation 
of these properties has been proposed, based on a characterization of 
high-noise and low-noise orderings. Furthermore, the dependence of 
noise on transfer function parameters and scaling has been investi- 
gated. These results point to an algorithm for ordering cascade FIR 
filters which has been implemented and tested for filters with a wide 
range of values of transfer function parameters. In every case, the 
algorithm gave results within expectation which were deduced to be 
very close to the optimum. Justification for the success of the algo- 
rithm has also been given. 
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Slope Overload Noise in Linear Delta 
Modulators With Gaussian Inputs 


By L. J. GREENSTEIN 
(Manuscript received September 22, 1972) 


This paper derives a slope overload noise power formula for linear 
delta modulators having ideal integrators and Gaussian random inputs. 
Although the same problem has been treated by others, the present result 
is the only one applicable to all slope-following capacities and input spectra. 

Despite its singleness of purpose, the paper divides logically into two 
parts. In Part 1, a common element in all previously published results 
as used to derive a new slope overload noise power formula. This derivation 
is analytically rigorous and provides some useful insights, but pertains 
to a particular kind of spectrum and so is incomplete. 

The more universal result we seek ts derived in Part 2. The approach 
here is far less rigorous and amounts to approximating the influences of 
other kinds of spectra by modifying the result of Part 1. The final expres- 
ston contains four spectrum-related coefficients, for which simple formulas 
are given, and has an estimated accuracy of 1 dB for all cases of practical 
interest. Computed results are given for two important families of spectra 
and comparisons are made with previously published results. 


Part 1 


I. INTRODUCTION 
1.1 Objective 


This paper presents some new theoretical results on slope overload 
noise in linear delta modulators. In particular, we derive the mean 
power of the (unfiltered) slope overload noise at the demodulator 
output when the modulator input is a stationary Gaussian random 
process and the modulator feedback path contains an ideal integrator, 
Fig. 1. This same problem has already been treated by other investi- 
gators, notably Zetterberg,! Rice (with O’Neal),? Protonotarios,? and 
Abate, but several factors have prompted a reexamination. One is 
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Fig. 1—Linear AM codec with perfect feedback integrator. (a) Equivalent block 
diagram. (b) Input and feedback signals. (c) Model used. 


that some important disparities among the results of these separate 
studies remain unresolved ; another is that the existing formulas do not, 
either individually or collectively, pertain to all slope-following 
capacities and input spectra; and finally, few explicit clues are avail- 
able as to just how accurate each formula is and where it loses validity. 
In response to this situation, we have endeavored to find an expres- 
sion for slope overload noise power that is accurate for all slope- 
following capacities and input spectra of possible interest. The result 
reported here satisfies that objective. 


1.2 Nowse Descriptions and Definitions 


The idealized AM codec (coder/decoder) shown in Fig. 1 ex- 
emplifies the process we want to analyze: Every 7 seconds the in- 
put signal (2(t)) is compared with a locally quantized version of 
itself (y(#)), and a unit impulse is generated with a polarity that is 
Gee if € 7 a The resulting binary impulse stream is 
applied through a feedback gain factor (6) and ideal integrator to 
produce y(t). At the decoder, y(é) is reconstructed by using a replica 
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of the coder feedback network, and a final low-pass filter smooths the 
sharp edges of y(t), yielding a closer approximation to a(t). 

Since y(t) cannot change by more than 6 units in 7 seconds, 6/7 is 
the highest input signal rate-of-change that the AM codec can follow. 
We call 6/7 the AM slope-following capacity and denote it by 73. When 
|dx/dt| exceeds this quantity, slope overload occurs and gives rise to 
the kind of error shown in Fig. 1b. In addition to this sporadic form 
of distortion, the “‘hunting”’ of y(t) for z(t) by means of quantum steps 
gives rise to a perpetual distortion called granular noise. Obviously, 
granular noise is reduced by decreasing 6, but at the expense of a 
reduced slope-following capacity and, hence, greater slope overload 
noise. 

To delineate slope overload and granular noise for purposes of this 
analysis, let us suppose that 2(t) is passed first through a AM codec 
having infinitesimal 6 and 7, but with the ratio between them the same 
as in the actual system (Fig. la). The decoded output (before the 
final filter) will then be the smooth function ¥(¢t) shown in Fig. le. 
If ¥() is then passed through the actual system, the decoded signal 
will be a very close approximation to y(t), Fig. 1b. In agreement with 
previous conventions, we define slope overload noise to be the differ- 
ence between x (é) and ¥(t) ; and the remaining distortion, (¥(t) — y(d)), 
is defined to be granular noise. This essentially is the approach tacitly 
followed in most of the published literature on AM noise.!~7 

We define a slope overload noise burst to be the nonzero difference 
between x(t) and ¥(¢) over an interval [t,,t)], which starts because 
|dx/dt| > 2, at t = t. and ends because ¥(t) intersects x(t) at t = b,, 
(e.g., see Fig. 2, in which ¢, = 0). The mean power of these bursts, 
averaged over all time, is the slope overload noise power. 

Two qualitatively different kinds of noise bursts can be identified. 
One is the kind initiated when |d2z/dé| increases through 2j while 7(t) 
is following x(t), which we call primary noise (see Figs. le and 4b). 
The other arises when a prior burst terminates at a point where | dx/dt| 
already exceeds xj, (see Fig. 4c) ; the new burst that commences at this 
point we call secondary noise. 

Finally, the slope overload factor is defined to be the ratio of the 
slope-following capacity to the rms input signal derivative, 1.e., 


A cos 


yr ma ee 
(dx/dt) 


(1) 


It should be obvious that the noise power decreases monotonically 
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Fig. 2—Zetterberg’s noise model. 


with increasing S. It is also clear that when S is large (in which case 
|dx/di| rarely exceeds x/), secondary bursts are rare so that primary 
noise dominates the noise process; and that, by the same token, 
secondary noise dominates the noise process when S is small. 


1.38 Equivalent Input Process 


The generality of our analysis will be enhanced if we can assume 
that {2(t)} is a bandlimited process. This assumption is made valid 
by regarding {x(t)} as an equivalent process related to the true one as 
follows: Let the true input process be {x,(¢)}, having a power spectrum 
X o(f). Since the delta modulator really acts on discrete samples of 
the input separated by 7 seconds, its response is the same as if the 
input were a process {a(t)} having sample functions of the form 


< sin (r(é — kr)/7) 
i) = okt) ——____— 2 
o) = Falke) (2) 
and a power spectrum 
2 1 
xin= Ex(s+"), Osfssony. — @) 
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We will assume here that the sample function and slope overload noise 
depicted in Fig. le correspond to the bandlimited process {2x(t)}, i.e., 
that the input spectrum is X(f). If {x.(t)} is bandlimited by some 
frequency W S 1/(27), the two processes, and consequently their 
spectra, are one and the same. (This condition is usually tacitly— 
though not always rightly-assumed in AM noise analyses.) If {2.,(t)} 
is not so limited, the difference between the two processes is aliasing 
distortion, which can be analyzed separately. 


1.4 The Spectral Moments 


The spectral moments of the input process play a decisive role in 
establishing the past and present slope overload noise formulas. 
Following the convention of most authors, we define the nth moment 
to be 


— i w?X (fds. (4) 


The bandlimited nature of X(f), as discussed above, guarantees the 
finiteness of all the moments. It is easy to show that the complete 
set of these finite moments, (bo, bi, bs, ---, bn, -+°), uniquely deter- 
mines X(f), just as the converse is true. 

It should be obvious that b, is the mean power of the nth derivative 
of x(t), i.e., bn = (x (2))2. Thus, for example, Vb; is the rms derivative 
of the input and so (1) can be rewritten as 


S = 2,/ Voy. (5) 


Also, bo is nothing other than the input signal power, as distinct from 
the ac input power which we denote by o?. If the input contains no 
discrete de component, then by = a?; otherwise, bo exceeds o? by the 
amount of the de power. 

Finally, we note from (4) that 6, is real and non-negative because 
X(f) is. An additional relationship for 6,, derived by applying the 
Schwarz inequality to (4), is 

bat 


b,2 ‘ n 
Dn—2 





IV 


2, (6) 


that is, b, for n 2 2 has a lower bound determined by the previous 


two moment values. We shall utilize this relationship later. 
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1.5 Outline of the Work 


Our objective is to find a suitable approximation to the true slope 
overload noise power formula, the latter being denoted by N(S). 
Section II describes the methods and results of previously published 
analyses. While the various approaches differ, they all suggest that 
this noise power is essentially determined by S and the first three 
spectral moments, bo, bi, and be. Section III treats this possibility as a 
premise and identifies a simple two-band process, the parameters of 
which can be chosen to yield any (bo, 61, b2) and which can be analyzed 
precisely. The slope overload noise power for this process is derived 
and is denoted by N,:(S). 

Unfortunately, N.,(S) is not precisely applicable to all spectra 
having the same bo, bi, and be. Instead, it is a lower-bound variation 
for all such spectra, which becomes increasingly unreliable in general 
as S decreases. To obtain a more universal estimate [denoted by 
N .(S) ] it is necessary to determine how much the slope overload noise 
power deviates from N,,(S) for spectra other than the two-band kind. 

To this end, Section IV derives a noise power formula, N.(S), 
applicable to all spectra but confined to the large-S region (S = 3.5), 
and Section V derives a noise power formula, N (8), applicable to all 
spectra but confined to the small-S region (S ~ 0). Section VI then 
combines these results to obtain a two-region approximation, N ,(S), 
which is accurate for all X(f) and S. The result is given by (62) and 
(77), where a, a2, Q1, Q2, a3, and a, are related to the spectrum parame- 
ters bo, bi, be, bs, bs, bs, X(O), and S? w2[X(f) — X(0) ]df. Further 
study shows that (77) alone can be used over all S with an accuracy 
of 1 dB or better up to S = 6.5, beyond which point N(S) is at least 
119 dB below the input signal power. For practical purposes, therefore, 
our final approximation to N(S) is just (77), with ai, de, a3, and a, 
given by (72) through (75). 

Section VII demonstrates the new result for two important families 
of spectra, and compares N,(S) with the noise power formula of 
Protonotarios. The latter is found to be highly accurate for most 
spectra of practical interest and S 2 2.0. 


II. REVIEW OF PREVIOUS WORK 


In the approaches of Zetterberg, Rice, and Protonotarios, N(S) is 
approximated as the product of the mean energy (&) per slope over- 
load burst and the average rate of occurrence (R) of such bursts. A 
burst is assumed to commence whenever |2’(¢)| increases through 
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the value xj, which means that only primary noise is considered in 
these studies. The average rate of these events is known from the 
earlier work of Rice® to be 


R= 1 esp { — (a)?/2b1}. (7) 


To find &, Zetterberg uses the model shown in Fig. 2, where the solid 
line of slope x4, represents the demodulator output for the duration of 
the burst, i.e., from ¢ = 0 to t = ty. Clearly, the noise is 


n(t) = x(t) — [a(0) + 2%t]; 0<t St, only (8) 


which is the first burst depicted in Fig. 2. To avoid deriving the 
random quantity ty, Zetterberg regards the slope overload burst to be 
the entire excursion of 2(t) above the line [2(0) + 2%], which is the 
variation m(t) in Fig. 2 rather than just n(é). This is obviously an 
approximation, since m(t) can contain one or more spurious bursts 
that are not really part of the original noise burst, as seen. Such 
bursts, however, are low in both energy and probability of occurrence 
when S > 1. 

Zetterberg proceeds by finding, at each t, the average of m? over 
the noise burst ensemble, and then integrating over all time and 
multiplying by R, (7). If done correctly, this leads to an estimate of 
N(S) which is marred only by the inclusion of spurious burst con- 
tributions and the ignoring of secondary noise effects. In fact, however, 
Zetterberg errs in defining the initial conditions of the burst and in 
averaging over these conditions, leading to an incorrect solution.* 
In addition, Zetterberg makes a number of functional approximations 
and at least one algebraic error.* For all these reasons, his final result 
is incorrect and will not be repeated here. 

The approach of Rice reverses that of Zetterberg, in that the energy 
per burst is found first and then averaged over the noise burst ensemble. 
The problem of spurious bursts is thus avoided, but at the expense of 
having to find t,. Rice accomplishes this by expanding x(t) into a power 
series about ¢ = 0, and then ignoring fourth-order powers in ¢ and 
higher, an approximation that gains in validity with increasing S. 
In addition, he approximates the third derivative of x(f) at t = 0 by 
its conditional mean given that 2’(0) = x/. This can be shown to be 
— (b2/bo)x4, and the relative variation of x’’’(0) about this value tends 
to be small when S > 1. 


* See Protonotarios® for a discussion of these errors. 
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Using these two approximations, and correctly identifying and 
averaging over the initial conditions of the burst, Rice obtains the 
result 


243 /bi\ 1 
Na(S) = (=*) = exp (—S?/2). (9) 


Note that the ratio of noise power to signal power (bo) depends solely 
on S and y £ bj/bobs. From (6) we see that y S 1 and observe that 
y = 1 only for an infinitely narrowband spectrum. 

It should also be noted that, as S—0, Nr(S) increases without 
bound at a rate S~—*. This is clearly not a true representation since 
N(S) should approach the ac signal power (oc?) in the limit of zero 
slope-following capacity.* 

Noting the disparity between Zetterberg’s result and Rice’s, even 
at large S where both should converge to exactness, Protonotarios has 
attempted to settle the issue by combining Zetterberg’s more accurate 
model with Rice’s correct averaging procedure. His analysis leads to a 
double integral solution which is exact except for the inclusion of 
spurious contributions and the ignoring of secondary noise effects. 
In reducing this formal result, Protonotarios makes some functional 
approximations and obtains the following: 


Nie () sen (=88/2VAGD (10) 
EN Aloe Nig) Ga . 
where 
25/8 2 1/8 
x = S/(bi/baby) (11) 


and A(-) is a function involving powers and exponentials of the 
argument, [eq. (66) of Ref. 3]. Once more, the ratio of noise power 
to signal power depends solely on S and y. In the limit as S > ~, A (x) 
approaches 1 and so N p(S) converges to Rice’s result, Nr(S). In the 
limit as S — 0, A(x) varies as S‘ so that Np(S) varies as S~ rather 
than S—5. This is still an unbounded increase, however, so that (10) 
is still not acceptable as a complete characteristic. 

Abate’s derivation of N(S) follows an approach quite different from 
the others. Using simulation results reported by O’Neal? for three 
particular spectra, he has developed an empirical relationship between 
the AM sampling frequency and the slope overload factor for which 


* If the input contains a discrete dc component, the delta modulator will follow 
it exactly so long as S ¥ 0; hence the result limso N(S) = o. 
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granular-plus-slope overload noise power is minimized. Combining 
that relationship with a simplifying approximation to van de Weg’s 
formula for granular noise power,® Abate obtains the following expres- 
sion for slope overload noise power: 


8r? by 

27 (27W)? 
where W is the signal truncation bandwidth.* We see that N4(S) goes 
to a finite value as S > 0, and that its variation at large S is expo- 


nential; in both respects, it differs from Nr(S) and N p(S). If we now 
define 


Na(S) = (1 + 38S) exp (—38) (12) 


A 8r? bo 


~ 27 (2nW)%by' we 





(12) reduces to the form 
2 
b 
Na(S) = K (~) (1 + 38) exp (—858). (14) 
2 


Evaluating K for the spectra studied by O’Neal, we find that it lies 
between 1.074 and 1.75 for the three cases, a range of just 2 dB. 
It is tempting to speculate, therefore, that (14) is a more correct form 
in general than (12), with K a universal constant on the order of unity, 
and that the apparent 2-dB spread in K over the three spectra is a 
result of experimental uncertainties. If we accept this notion, we again 
have the result that the ratio of noise power to signal power depends 
solely on S and y. 


III. NOISE POWER IN TERMS OF Do, by, bo 


3.1 The Two-Band Process 


The published results cited above suggest that N(S) is essentially 
the same for all processes having the same values for the zeroth, first, 
and second moments. Assuming for the present that this supposition 
is correct, we call attention to the two-band process, the spectrum of 
which is shown in Fig. 3a. In the limit as the two bands become in- 
finitely narrow, the zeroth, first, and second moments become pre- 
cisely bo, bi, and bs. By analyzing this simple process, therefore, we 
can derive an exact noise power formula [N.s(S) ] in terms of bo, bi, 
and by and then see how universal it really is. 
 *The three spectra treated by O’Neal are of the truncated Butterworth type, 


with corner frequency-to-bandwidth ratios of 0.068, 0.25, and «. The data used to 
derive (12) cover an S-range from about 1.4 to 4.2. 
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x (t) 





Fig. 3—The two-band random process. (a) Power spectrum. (b) Sample function. 


To appreciate the simplicity of the two-band process for purposes 
of analysis, we should view it in the time domain. A sample function 
is shown in Fig. 3b and is seen to consist of a dc level, D, plus a sinusoid 
having radian frequency w,. = Vb2/bi, phase ¢ = wots, and amplitude 
A. Dis a Gaussian variate, whose mean and variance across the sample 
function ensemble are 0 and (bo — b7/b2), respectively; ¢@ is a uni- 
formly distributed variate on [—7z, +7]; A is a Rayleigh variate of 
mean-square value 2b7/b.; and D, ¢, and A are mutually independent. 

We can derive the slope overload noise power for this process by 
finding the mean noise power associated with a given sample function 
and averaging over the distributions on D, ¢, and A. This approach is 
simplified by the fact that D and ¢ do not really influence the result. 
That is, the noise pattern for a given sample function converges 
ultimately to a variation about D that depends solely on the sinewave 
amplitude and frequency (Fig. 4). 

We thus see that the noise power per sample function is the steady- 
state noise energy per half-cycle, denoted by &(A), times the rate of 
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_ x(t) AND y(t) 


~- x(t) AND y(t) 





Fig. 4—The three regions of sinewave amplitude. (a) 0 S$ A S 24/w. (no noise). 
(b) (x}/wo) <A < (xh/w.) V1 + 72/4 (primary noise). (c) A 2 (x'/wo) V1 + 72/4 
(secondary noise). 





occurrence of half-cycles, which is Vb2/b1/n. The total noise power for 
the two-band process is then the average of this quantity over the 
distribution on A, i.e., 


e ~ i yp 
ae by 


0 


&(A)p(A)dA. (15) 
SY” 
(pdf of A) 
Since A is known to be a Rayleigh variate of mean-square value 


(2b7/b.), the only unknown is the energy function 6(A). 


3.2 Derwation of &(A) 
Three distinct regions of A can be identified, each giving rise to a 
distinct pattern of signal and noise. The region depicted in Fig. 4a 
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wy 3 &/ (x'g)” 





—1.0 -0.5 0 0.5 1.0 1.5 2.0 qr2 3.0 3.5 
€=(Awo/ x‘)? ~1 


Fig. 5—Noise energy per half-cycle as a function of A. 


(Region 0) is one in which no slope overload occurs because |2’ (t) | 
< 2} at all t. Clearly, 6(A) = 0 in this case. The region depicted in 
Vig. 4b (Region 1) is one in which slope overload occurs in each half- 
cycle, but over less than the complete interval. We see that the noise 
in this case is primary noise, as defined earlier. Finally, the region 
depicted in Fig. 4c (Region 2) is one in which slope overload occurs 
for the entire duration of each half-cycle. The noise in this case is 
seen to be secondary noise.* 

To find 6 in Region 1, we analyze the burst spanning [1#,, t. + At] 
in Fig. 4b. The noise energy in this burst is 


totAt 
6 = | [A cos wot — {xo — xo(t — t.)} ]?dé. (16) 
to 


The quantities 2, and ¢, are found by means of trigonometric identities 
to be 





1 | X', 
to = VA? — (25/w,)? and ¢,= — sine ( ). (17) 


Wo Wo 


Unfortunately, the quantity At cannot be solved for explicitly, but is 


*For the special process under discussion, we see that a given sample function 
has either no noise at all, primary noise only, or secondary noise only. For more 
spectrally distributed processes, both kinds of noise occur in the same sample func- 
tion and in a less regular pattern. 
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related to A by the transcendental equation 


w At — sin (w,At) 7? Aw,\? A 
fe - ( ) ee (18) 
1 — cos (w,At) 3 


We can relate & to A (or £) by combining (16), (17), and (18) and 
eliminating the common parameter, w,At, between (16) and (18). 
The result is a unique correspondence between w3&/(x/)? and & having 
the variation shown in Fig. 5 for Region 1. (Note that Region 1 cor- 
responds to the &interval [0, 7?/4].) A highly accurate functional 
description of this variation has been found to be 











3 
Wo8 81 
= £7/270.43 exp (—0.5014£) + 0.57 exp (—2.108) ]; 
(x)? 1407 
OS E< 77/4. (19) 


Comparison with exact results shows this function to be accurate to 
within 0.8 percent. 

To find & in Region 2, we analyze the burst spanning [¢,, t. + At] 
in Fig. 4c. In this case, x, and t, are given by (17), as before, but 
w, At is precisely a for all A. Applying these relationships to (16) and 
performing the indicated integration, we obtain 


cpral-0-a)} ees 


This variation is also shown in Fig. 5. 





3.3 Expression for N1»(S) 


We now have expressions for & in the three regions of A and can 
average this complete result over the distribution on A (or &) to find 
Ni», (15). By using (18), along with w? = b2/b,, S = 25/Vbi, and the 
fact that A is a Rayleigh variate of mean-square value (2b{/b:), we 


obtain 
S? S? S? 
p(é) = woe (- =) exp (- = 8); é=-1 


= 0; elsewhere. (21) 
Combining this with (15), (19), and (20), we obtain 
bi bi bi 

NalS) = (=) Fy(S) + (~*) ris) =(-*) (8) 22) 
be be be 


Primary Secondary 
Noise Noise 
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where F'(S) & Fi(S) + F.(S), and 


F,(8) = ~— (=) s exp (- =o. 43Q = (> +0. 5014)| 


+ 0.57 = (> ake 2.10)} | (2) 





F,(8) = E 4 (~—)s*| exp | - (1 - =) >: (24) 
24 4/ 2 
Q{u} = aa [ron et (vu) 
ysl? 16 


7 35 105 
= jure as a nie = a8 = wn] exp (-w)| > “25) 


The variation F (S) is shown in I’ig. 6 along with its component parts, 
F,(S) and F.,(S). These curves indicate, for the two-band process, 
which region of S is dominated by primary noise and which region 
by secondary noise. 

The only inexactness in our result for N,(S) arises from the func- 
tional fit to 6 in Region 1, (19). Since this fit is accurate to within 
+0.8 percent at all &, and both &(£) and p(&) are non-negative func- 


\ 
F,(S), (PRIMARY \ 
NOISE) 


(F(S), (TOTAL NOISE) 


NOISE POWER IN dB ABOVE b? /b, 





SLOPE OVERLOAD FACTOR, S 


Fig. 6—Normalized noise powers for the two-band process. 
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NOISE POWER IN dB ABOVE (b? /bo) 


re (bt by by) = 0.11 


PROTONOTARIOS:: 
—- —(b? / by by) = 1.00 





SLOPE OVERLOAD FACTOR, S 


Fig. 7—Comparison of Niy(S) with previous results. 


tions, the given expression for N,,(S) is also accurate to within +0.8 
percent (or +0.035 dB). 


3.4 Comparisons and Interpretations 


The variation of VN,» with S is plotted in Fig. 7, along with the varia- 
tions of Nr, (9); Np, (10); and Nu, [(14) with K = 1]. Of particular 
interest is the fact that N.,(S) converges with the results of Rice and 
Protonotarios as S— «. The approximations underlying the analyses 
of Rice and Protonotarios become increasingly valid as S increases, 
so the convergence of N..(S) with their results is not surprising. The 
apparent dependence of noise power on the moments b; and b2 alone 
at large S is explained by the fact that slope overload bursts are of 
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short duration in this region; hence, they are shaped primarily by the 
lower-order curvature of z(t), which is reflected in b, and bo. As S 
decreases towards small values, however, the higher-order curvature 
of x(t), reflected in b;, bs, etc., also influences the noise bursts and, 
consequently, the noise power. We should therefore expect the true 
noise power [N (8) ] associated with a given spectrum to reflect other 
features of that spectrum at low S besides b; and b». If this is so, then 
N.»(S) cannot be assumed to be a universal formula. 

To show that this expectation is correct, we note that, as S — 0, 
N.» converges to b{/bs, which is precisely the ac power (2) in the two- 
band process. This is a quite general result, i.e., N(S) always con- 
verges to ao? as S > 0. However, some other process having the same 
bo, bi, and b, can have an ac power anywhere between b7/b. and bo, 
so that N(S) will not converge to bj/be for all spectra having these 
moments in common. 

We conclude, then, that N(S) and N.,(S) converge at high S but 
become increasingly dissimilar as S decreases towards zero, the dis- 
similarity depending on the differences between the actual process 
spectrum and the two-band spectrum having the same bo, bi, and be. 
In Part 2 (following), we will derive a factor which relates N,(S) 
and N(S) when X(f) is not a two-band spectrum. 


Part 2 
IV. NOISE POWER FOR S = 3.5 


We now derive a noise power formula applicable to all spectra, with 
S large (primary noise only). The derivation combines the power series 
method of Rice with the infinite-time averaging procedure used by 
Zetterberg and Protonotarios. The differences here are that many 
terms of the power series expansion are included, and steps are taken 
to minimize the contributions from spurious bursts. The complicated 
expression that results is reduced to a simple formula [N1(S) ] that 
displays explicitly the influence of the higher-order spectral moments. 


4.1 Power Series Representation of the Noise Burst 

Let t = 0 be the time origin of a particular positive-going burst, so 
that x’(0) = x, and x’’(0) > 0, e.g., Fig. 2. Beginning at t = 0, ¥(é) 
follows the straight line «(0) + a(t until this ramp function intersects 
x(t) again. We can therefore define the quantity 
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v(t) 2 x(t) — [e(0) + xt]; 2 >0 (26) 
and write the variation of the noise burst as 
n(t) = v(t) ifv@’) >0 forall O<t St 
= 0 otherwise. (27) 


Henceforth, we shall use the symbols n; and v; to denote n(t) and v(é), 
respectively. 

Following Rice, we can express x(t) by a power series expansion 
about ¢ = 0. Noting that x’(0) = xj and denoting 2’’(0) by 2)’, v; can 
be written as 


1 ] 
Ve = a0 + — a" (0) + --- —a™(O)i™ + ---, (28) 
3! m! 


We now invoke the property of Gaussian random processes that the 
mth and (m — 2)th derivatives at a given time instant are related by 


m—1 





rim) = — ged, (29) 


m—2 


where dm is a zero-mean Gaussian variate of mean-square value 


“2 Dm-1 


es eee 





(30) 


With this relationship we can write v, as follows: 


| E bs t4 bs bs t® | 
Se as Se ee ee 
2! by 4! be by 6! 
E f2 ba by a dy be EI | 
— Lo Cece i ee eet 
bi 3! bi bs 5! by bz bs 7! 


+ {2 ea Ee a]? E a)o+ | 
ig aki? 9 om lo ae ee ema 1) 


Ait 6 (31) 


where A, is the first bracketed quantity and 6, is the second. We see 
that A, is the conditional mean of v, given that 2x'(0) = 2 and 
x’ (0) = x4’; and that 6, is a linear sum of mutually independent zero- 
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mean Gaussian variates (d3, ds, etc.) whose mean-square value at 
time t is therefore 


“2 2 t by 08 a8 
i: ([-2a+- Jat] op 
4! 0b, 6! 
We thus have the result that the conditional pdf of v at time ¢ is 








p{o.|x%, ot = 


1 1 v.— A, 2 
ar? [5 . yf 2 


4.2 Derivation of Mean Noise Power 


Let n? be the conditional mean square value of n at time ¢, given 
that x’(0) = x3 and x’’(0) = 2// > 0. To find the mean burst energy 
when x’(0) = 2%, we must average n? with respect to 2/’ and integrate 
the result over all time. We can then multiply this mean energy by the 
mean rate of primary noise burst occurrences to obtain the mean noise 
power. This approach is valid so long as secondary noise can be 
neglected, which it can in the region S = 3.5 under consideration. 

From Protonotarios,? we know that the correct density function for 
averaging with respect to the conditions x’(0) = x3 and xi’ > 0 is 





ney = x! ae |- (2’}) r al! > 0 (34) 
be Dba)” _— 


while from Rice® we know that the mean rate of noise burst occur- 
rences is (Vb2/bi1/7) exp (— S?/2). The noise power in the region S = 
can thus be accurately given by 


n9-Eo(-2)f [2 


| (xo)? 
-exp |— 


2be 





n?(x', x0 )dxvdt. (35) 


The remaining analytical task is to find n; n2 (al, 08 ah 

The method of Protonotarios in this regard amounts to averaging v7 
over positive values using the conditional pdf of v,, (33). To do this, 
however, is to incur the spurious burst problem depicted by Fig. 8a. 
The following line of reasoning leads to an improved procedure. 
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Fig. 8—An illustration of v, and n:. (a) 1 and n:. (b) vz and v}*. 


Let A{” denote the partial power series for A; up to and including 
the fourth-order term, and let 


po A ee (36) 


We can show that v{” =~ v, when tis small and that v, > v{” in general. 


What we have found in particular is that, for S greater than about 2, 
the difference between v, and »;* does not widen appreciably over a 
typical burst duration (see Fig. 8b). Also, since A;” is a strongly de- 
creasing function of ¢ beyond its peak, v{* tends to be the same. 

We therefore calculate nj according to the following criterion: A 
given value of v; is counted as part of the noise (i.e., n; is assumed to 
be v,) if either {v, > 0, vf > 0} or {v, > 0, vf < 0, % < 0} ;* other- 
wise, v; 1s not counted and n, is assumed to be 0. For the sample func- 
tion shown in Fig. 8b, this means that the energies in intervals A and 
B are both counted in the calculation of n?; the energy in interval C 
is not counted, since v{ < 0 and d, > 0; and the energy in interval D 
is counted, even though it is clearly spurious. We can reduce contribu- 


* The symbol #; stands for dv(t) /dt. 
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tions of the latter type, however, by using a suitably finite upper limit 
in the integration over time in (35). That is, we choose an upper limit 
sufficiently high that virtually all legitimate noise contributions are 
counted in the calculation but sufficiently low that some spurious 
contributions are omitted. By studying the time variations of A, and 
6, under a wide variety of conditions, we have settled on an upper 
time limit of 3.5 Vby/b. 
Using the above criterion, we have been able to show that 


n? = > (1 + 2a?)[1 + erf (a4) ] +e (2a — a4) exp (-a)} 
+ = (erfe (4)) }(1 + 2a2)[erf (a) ~ erf (a:)] 
+ = [a exp (—a?) — (2a — as) exp (—a)}] (37) 
where: 
a = A,/ Wah; (38) 
an = At /1BaF: (39) 
and 


a 2 A,/f2(5)2. (40) 


To complete the derivation of N(S), we insert this complicated expres- 
sion into (35) and perform the integrations over x,’ and t, remembering 
to use an upper limit of 3.5 Vbi/bs in the time integral. 


4.3 Simplified Formula 


The formal solution just described is believed to be more exact 
than that of Protonotarios [eq. (53) of Ref. 3], because it entails a 
substantial reduction in the spurious burst contributions. Furthermore, 
in reducing his formal solution to computation, Protonotarios makes 
a number of functional approximations which obscure the influence 
of the higher-order moments 63, b4, etc. In reducing the present solu- 
tion to computation, we must also make approximations but of a 
different kind, i.e., those associated with performing numerical double 
integration and truncating the infinite power series in the integrand. 
The impact of the higher-order moments, however, is not approxi- 
mated away in the process; in fact, the following discussion develops 
an empirical expression [denoted by N,(S)] which exhibits these 
influences explicitly. 
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We define a dimensionless quantity, which we call the nth-moment 
exCess, as 


A bn — (Dn—1/Bn—2) 


M,, : eS (41) 
(B9-1/5"-2) 


From (6) we know that 17, 2 0 for all n 2 3 and, using Fig. 3, we 
can show that MV, = 0 for all n 2 3 if and only if X(f) is the two- 
band spectrum. We can therefore say that N(S) differs from N..(S) 
only to the extent that M3, M,, etc., are not zero. A simple model 
that reflects this fact is one which describes the logarithmic difference 


between N and N;, as a linear combination of the moment excesses, 
In(N/N i») = AM3;+ BM,+CM;+4+ --: (42) 


where the coefficients A, B, etc., are, in general, functions of S. To 
the extent that this representation is accurate, it is reasonable to 
assume further that the higher-order terms, involving M5, Mz, etc., 
are negligible for S = 3.5. The reason is that these moment excesses 
have little influence on slope overload bursts at such high slope- 
following capacities. The result of this line of reasoning is an approxl- 
mation of the form 


N(S 2 3.5) + Nis(S) exp{A(S)Ms + BUS)Ms + C(S) Ms} 
£Ni(8). (48) 


We have tested this approximation by applying the double integral 
solution for N(S), (35) and (87), to a number of spectra, using a 
ninth-order power series to represent v, With the computed noise 
powers for these spectra, we have then obtained results for A(S), 
B(S), and C(S) over the range 3.5 S S S 10.0 (Fig. 9) which appear 
to be quite accurate. In particular, the use of (48) with these results 
is found to predict N(S = 3.5), as computed from (85), to within 
0.2 dB for all spectra of practical interest. 

There exist, however, special conditions on X(f) for which (43) 
yields too high estimates of the true noise power. This can occur when 
X (f) contains an isolated high-frequency component which contributes 
materially to 63, bs, and/or bs but not to bo, b;, or be. Under such circum- 
stances, x(t) might change too rapidly during the noise bursts to be 
properly analyzed by the power series approach as used here. Calcu- 
lations indicate that these cases fall outside the range of practical 
interest, so we need not consider them further. 
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Fig. 9—Empirical results for A, B, and C; S 2 3.5. 


V. NOISE POWER FOR S =~ 0 


We now derive another noise power formula applicable to all spectra, 
but with S very small (secondary noise only). The representation we 
seek is a partial power series of the form 


N(S = 0) = o[Co + CS + C.S2)] & N,(S) (44) 
where o” is the ac power of the input process. 


5.1 Formulation 


The analysis assumes {2x(t)} to be just the ac part of the input 





process, so that x?(t) = o?. Since de components have no effect on noise 
power for S ¥ 0, we are justified in ignoring them. 

We consider x/ to be so small compared to Vb, (i.e., S so small 
compared to unity) that the delta modulator is always in slope over- 
load. In this case, the decoded output [denoted here by y(t) rather 
than ¥(é) ] is just a succession of alternating ramps of slope +24, the 
slope polarity reversing whenever y(t) intersects x(t). This situation is 
illustrated in Fig. 10b. A valid model of the delta modulator for 
analyzing this case is given by lig. 10a (ef. Ref. 9), from which we 


see that the noise signal is 
t 


n(é) = x(t) — y() = 2(t) — 2% | POY (45) 


—o 


where 
nq(-) & sgn{n(-)}. (46) 
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x(t), y(t) 





Fig. 10—Equivalent delta modulator and signals for S ~ 0. (a) AM model for 
S = 0. (b) Input and feedback signals. 


The noise power in the region S ~ 0 can thus be accurately given by 
N(S) = x(t) — 22’, Eo [ nal) 
—2x(d)y() 
+ (x)? [ i: nq(u)n_(v)dudv. (47) 
~o Jw 
y>(2) 


If we envision each of these three terms as a power series in x), we 





see that x?(t) has a zero-order term only (i.e., is independent of 2) ; 
—2zx(t)y(t) may contain first-, second-, and higher-order terms; and 





y?(¢) may contain second-order terms and higher. We can thus rewrite 
(47) as 


N(S) = Ko + [Kiva + Kea(xe)? + +++] + [Kas(x0)? +--+]. (48) 
~—— Fe ee 


x*(t) —22(t)y(t) y(t) 


Retaining only terms up to second-order, substituting 2} = Svb,, and 
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comparing with (44), we see that 
Ko K,Vby C, = (Koa + Kos)bi 


Co =—;} 1= ; 2 
o? o? 





— (49) 
Clearly, Ko = x(t) = o’, so that Cy = 1 in all cases. The remainder 
of this section is devoted to finding C,; and C2. 


5.2 Derivation of Cy 


We will derive C; by evaluating the first-order term of —2z(t)y(t), 
(47), and then invoking (49). For x; very small, n,(t) does not differ 
appreciably from x,(t) = sgn{2(t)}, which is the quantizer response to 
x(t) applied alone. It is therefore fruitful to represent n,(t) as 


Nq(t) = X(t) + A(t) (50) 


where 6(t) is a correction signal related to the finiteness of x) and is 
defined by 
+2 if y<x<0 
6=.+-2 if y>ux>od (51) 
0 otherwise. 


The relationships between n,(t), x,(t), and 6(¢) are illustrated in Fig. 
11. Although the peak magnitude of @(¢) is fixed, its “duty cycle” can 
be seen to depend directly on 2%. 

Inserting (50) into the cross term of (47) and interchanging the order 
of integration and averaging yields 


t 


—2a(Hy(t) = — zz | 


oO 


x(t)a(u)du — 22% [ x(t)@(u)du. (52) 


Since the duty cycle of 6(¢) vanishes as x, 0, we can see that the 
first term must be identical to K,x; in (48), while the second term 
contains K»,_(2z,)* plus higher-order terms in 2%. 

To evaluate the first term, we use a relationship applicable to 
{x(t)} because it is a Gaussian process, namely, 


where R,(r) is the autocorrelation function of {x(t)}. We now apply 
this to the first term of (52) and use the fact that the area under 
R.(r) from r = 0 to r = © is 1X(0). Equating the result to Kix, 
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Fig. 11—The relationships between x(t), y(t), x¢(t), ng(t), and A(t). 





we obtain 
gerne 1 X(0) (54) 
‘NOR og 
which, from (49), yields 
a 1 fx(0) a (55) 
oo V20 | a 


5.3 Derivation of C2 


The general approach used in the preceding analysis was also applied 
to finding C,. Unfortunately, it leads to an approximation which is 
both highly complicated and not sufficiently precise for many spectra. 
For this reason, we have resorted to finding C, for a not-quite-Gaussian 
process having the same spectrum as the specified Gaussian one, and 
verifying that it is virtually exact under what appear to be “worst- 
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case’ conditions. The details, which are documented but not pub- 
lished,!° are quite involved and so we shall merely outline the approach. 

Consider a wide sense stationary process {X(t)} having sample func- 
tions of the form 


X(t) = A cos(at + ¢) (56) 


where A is a Rayleigh-distributed amplitude of mean-square value 
2c?; ¢ is a random phase uniformly distributed over [—7z, +7]; and 
@ is a random frequency whose pdf is 


] 
p(@) = sa = ||/2m); alla. (57) 


It can be shown that the power spectrum for this process is X (f), and 
that the ensemble pdf of ¥ at any ¢ is Gaussian with zero mean and 
variance o?. Thus, {X(¢)} has certain properties in common with the 
Gaussian process {z(t)} having the spectrum X (f). 

We can derive the mean slope overload noise power for the process 
{x(t)} by finding the noise power for a single sample function, (56), 
and averaging over the distributions on A, ¢, and 6. Using the methods 
of Section IIJ, we can show that the averaging over A and ¢ yields 


Naia) = 1a" (—) (58) 

x, |@ 
VA% 

where F(-) is defined by (23) ) through (25). Averaging (58) over the 

pdf of a, (57), and replacing AZ by 20? and 2 by SVbi, we obtain 


bole 








: . iS 
N(S) = | x(pr ( ) ay. (59) 


2afo 

The derivatives of this function with respect to S can be found quite 

simply, although some caution is required when evaluating them at 
S = 0. The resulting value for C2(=4N’’(0)) is found to be 


comin (Ef aa Gee) | 


Tor most realistic processes, the continuous part of the power spec- 
trum has zero slope at f = 0, i.e, X’(0) = 0. In that circumstance, 
(60) reduces to the much simpler form 


Cy = (= -2)- ee 


- ome df; X0)=0. (61) 
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Fig. 12—Two families of spectra. (a) Truncated Butterworth spectrum. (b) Band- 
pass uniform spectrum. 


In general, the above C, for the non-Gaussian process {X(¢)} can- 
not be compared with the C, for the Gaussian process {2(t)} since we 
do not have a result for the latter case that is good for all X(f). We 
do, however, have exact results under some special conditions on X (f) 
for which C, is expected to be maximally dissimilar for the two pro- 
cesses. The very close agreement observed under these conditions 
persuades us that (60) is a reliable representation of C. for Gaussian 
processes having any X (f). 


5.4 Final Result for N (8S) 


Our final approximation to N(S ~ 0) is (44) with Co = 1 and C; 
and C, given by (55) and (60). Note that when X (0) = 0, C1 is zero 
and C; is negative, so that N,(S) is convex at the origin; and that when 
X (0) # 0, Ci is negative and C2 can be either positive or negative, 
depending on X(f). For a spectrum like the one in Fig. 12b, this 
implies a functional discontinuity at 82 = 0, which can be explained 
as follows: So long as 6, ¥ 0, the curvature of N,(S) is convex at the 
origin, but it becomes concave at some S > 0. As 2 decreases, the 
curvature at the origin becomes sharper and the point of inflection 
occurs at successively lower values of S. Finally, in the limiting case 
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B. = 0, the point of inflection occurs at S = 0 and so N,(S) becomes 
concave at the origin. Thus, an abrupt increase in X at f close to 0, 
reflected in a large negative value of C2, signifies a very small range 
of validity and calls for caution in using the new result for N(S = 0). 
For concreteness, suppose that the continuous part of X(f) in- 
creases to a large peak density at some low frequency f., where 
fo < Vbi/2xc. To obtain a more useful approximation to N(S ~ 0) 
under such circumstances, C'; and C, should be recomputed as if the 
point of high density were at f = 0 instead of f = f,, i.e., as if the 
spectrum were X(f + f.). The resulting NV.(S), (44), then corresponds 
to a spectrum quite similar to X(f) but without the abrupt change 
at low f. Consequently, it exhibits the correct ‘large-scale’? curvature 
at low S while omitting the sharp ‘‘small-scale’”’ curvature at S near 
0 associated with the true C; and C,. An empirically derived rule-of- 
thumb that leads to good results in all cases is the following: If 
C. < — 10, shift the continuous part of X(f) to the left until the 
nearest high-density peak is at f = 0, and recompute C; and C? for this 
modified spectrum. Then compute the quantity C = (0.5 C2 — 10 C1) 
using these values. If C’ exceeds the original value of C2, use the newly 
computed values of C; and C:; otherwise use the original ones. 


VI. GENERAL NOISE POWER FORMULA 


We now seek an expression for slope overload noise power [N.(S) | 
which is accurate over all S for all spectra and is simple to use. Our 
approach is to derive separate functional descriptions for the regions 
0<S <4.0andS 2 4.0. The descriptions are such that the resulting 
N .(S) and its first derivative are continuous at the boundary S = 4.0. 


6.1 NCS) for S 2 4.0 

The result Ni(S), (43), gives a very accurate approximation to N (S) 
in the region S = 3.5. However, our result for Ni.(S), (22) through 
(25), is quite complicated, and the results for A(S), BGS), and C(S) 
are available in graphical form only (Fig. 9). We now use these results 
to derive a simple but still accurate representation for this region. 

Taking a suggestion from Rice’s analysis,? we speculate that N (S) 
for S = 4.0 can be accurately approximated by 


243 bi exp (—S?/2) 
4V2rb, 88 
-Lexp (—a1 exp (—a28)) | e NS 2 4.0). (62) 


N(S = 4.0) = 
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The quantity preceding the brackets is Rice’s result Nr(S), (9), while 
the bracketed term (with az > 0) is a monotonically decreasing func- 
tion that converges to unity as S— ©. We choose a; and a: so that 


N (4.0) = Ni(4.0) and N.(4.0) = Ni(4.0). 


By using the data in Fig. 9 to obtain A(S), BGS), CGS), and their first 
derivatives at S = 4.0, and the data in Fig. 6 to obtain N..(S) and 
its first derivative at S = 4.0, we obtain 


0.2141 — 0.024M; — 0.1961, — 0.0671; 


= 63 
0.6394 — 0.057M; — 0.4261, — 0.093M; ne 


ae 
and 
a, = (0.6394 — 0.057173; — 0.426M, — 0.093M;5) exp (4a2) (64) 


where 3, M., M; are the moment excesses defined by (41). Note that 
when X(f) is the two-band spectrum, (1, = 0, n = 3), a: and a» 
reduce to 2.45 and 0.336, respectively. 

The above approximation to N(S = 4.0) can be tested by comparing 
it with the more precise results computed for various spectra using 
(85) and (37). The indications from such comparisons are that the 
approximation is accurate to within +0.3 dB over S = 4.0 for all 
spectra. 


6.2 N.(S) for0 < S < 4.0 


The analyses of Sections IV and V inform us about N(S) at the 
extremities of the region 0 < S < 4.0, but not about the variation in 
between. To estimate this variation reliably, it is convenient to express 
N(S) in the form 


N(S) = GS)N (8) (65) 


and then seek an approximation to G(S). The latter is a spectrum- 
dependent function that differs from unity because—and to the extent 
that—X (f) differs from a two-band spectrum having the moments bo, 
by, and by (Fig. 3). From physical reasoning given earlier, we expect 
that G(S) > 1 as S— o and that, as S decreases from high values 
towards 0, G(S) increases because—and to the extent that—b;, b., etc., 
exceed the minimum values, (6), corresponding to the two-band spec- 
trum. The maximum value of G(S) should then be its value as S — 0, 
which is o7(b2/b}). 

This reasoning is supported by careful scrutiny of the results for 
N.(S) and N,(S8), Sections IV and V. From Section IV we can show 
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that G(S) is a nonincreasing function* for S = 3.5 and that it goes to 
unity as S— o. From the results of Section V we can show that G(S) 
is a nonincreasing function* at S ~ 0 and that it goes to o7(b2/b%) as 
S— 0. A logical consequence of these observations is that, for0 < S 
< 3.5, G(S) is either a nonincreasing function, or has at least two 
extrema. The latter possibility has no physical basis in fact, and so we 
conclude that G(S) is a nonincreasing function over all S, approaching a 
maximum value of o°(b:/b?) as S— 0 and a minimum value of 1 as 
S— 0, 

Computations show that the decrease in G(S) over 0 < S < 4.0 is 
less than 10 dB for all spectra of possible practical interest. Thus if we 
can find a functional approximation to the quantity 


N(S) 
N.2(S) 





g(S) =n G(S) = in] i 0<8 <40 (66) 
that is accurate to within +10 percent, the resulting approximation 
to G(S) will be accurate to within +1 dB for all spectra of interest. 
This accuracy should be possible to achieve since g(S) near S = 0 
and S = 4.0 are known from the results of Sections III, IV, and V, 
and the above argument persuades us that g(S) is nonincreasing over 
the region in between. 
The function that we use to approximate g(S) is 


go(S) = In (o bo/b3) + ai + afexp (a8 +a) —1] (67) 


where a1, «+ a4 are chosen to give g,(0), go’(0), go(4), and g3(4) the 
values predicted by the results of Sections III, IV, and V. These values 
are found by using (66), with N,(S) and N.(S) replacing N(S) at 
S = Oand S = 4, respectively. Thus, the following equations must be 
satisfied : 

ay + aed3 = C1; (68) 


2 2 an 
a2(a3 + 204) = 2C2 — Cy + (4 = =); (69) 


In (o bo/b;) + 4a1 + aslexp (4a3 + 16a,) — 1] 
= 0.057M; + 0.426M,+0.093M;; (70) 
and 
a, + a2(a3 + 8a4) exp(4a3 + 16a4) 
= — 0.024M; — 0.196M, — 0.067M; (71) 


*In fact, G(S) is a decreasing function for all but the two-band spectrum, for 
which case it is precisely 1. 
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TABLE I—lF’oRMULA COEFFICIENTS FOR THE BUTTERWORTH SPECTRUM 





Bi ar a2 a3 a4 
0.02 —0.177 2.26 —3.93 —8.76 
0.068 —0.151 1.29 —3.70 —7.04 
0.10 —0.137 1.06 —3.75 —6.99 
0.25 —0.096 0.66 —3.90 —7.03 
0.50 — 0.066 0.50 —3.91 —6.67 
1.0 —0.047 0.42 —3.86 —6.22 
2.0 —0.040 0.38 — 3.84 —6.00 

eo — 0.036 0.37 —3.83 —5.90 














where Ci, C2. and M,(n = 3) are given by (55), (60), and (41), 
respectively. 
A set of solutions that is quite accurate in most cases is the following : 


a, = — 0.024M; — 0.196M, — 0.067M; (72) 
In (o bs/by) — 0.155M; — 1.210M,—0.361M; (73) 
(Cy — a1)/a2; a. > 0 
a3 = 0; a, = 0 (74) 
[C1 — ay| /a2; a, <0 
2 w 
2C2— C1 + (4- =) 
See eg te ee. «8 78) 


I 


Qe 


1 
a, = <- 
: 2 Qe 

0; a2 <0 

The results for a3 and a, when a2 S 0 are approximations only, but a» 
is nonpositive only when X(f) is relatively narrowband, i.e., when 
o°b,/b; ~ 1, in which case G(S) is close to 1 over all S. Hence, these 
approximations lead to a quite accurate representation of g(S). In the 


TABLE I[J—ForMULA COEFFICIENTS FOR THE UNIFORM SPECTRUM 





Be a1 Qe a3 a4 
0 — 0.036 0.37 —3.83 —5.90 
0.05 —0.036 0.32 —4.44 —8.21 
0.10 —0.036 0.27 0.14 —11.96 
0.20 —0.035 0.17 0.21 —7.58 
0.30 —0.033 0.08 0.41 —7.95 
0.40 —0.029 0.02 1.44 — 18.30 
0.50 — 0.024 —0.01 —1.76 0 
0.80 —0.004 —0.01 —0.43 0 
1.0 0 0 0 0 
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NOISE POWER IN dB ABOVE o2 


BUTTERWORTH SPECTRA 
——— UNIFORM SPECTRUM 





SLOPE OVERLOAD FACTOR, S 


Fig. 13—Noise power results for several spectra. 


more usual case where a2 > 0, the above solutions are highly accurate 
so long as 4a; + 16a, S$ — (4+ Inaz), a condition satisfied for most 
spectra. If this inequality is violated, solutions for a;, --- a4 can be 
found by graphical or computerized methods, but such situations 
appear to be rare. 

With g(S) approximated as above, N(S) over 0 < S < 4.0 can be 
expressed as the product exp (g.(S))Nis(S). To obtain a usable ex- 
pression, however, Ni.(S) should be given in a more simple form than 
eqs. (22) through (25). An approximation to N,.(S) which is accurate 
to within 0.3 dB is 

2 


b 
Nis(S) = = [1 + 2.7538 + 2.95287] 
2 
-exp (—2.7538S — 0.34182; O0<S<4.0. (76) 
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Fig. 14—Noise power results in greater detail. 


Combining (76) with (67) leads to the following: 


N..(S8) = o°[1 + 2.7538 + 2.95282] exp (—0.34182) 
-exp {(a1 — 2.753)S + a.[exp (a38 + aS?) — 1]}; 
0<S < 4.0. (77) 


6.3 Final Expression 


A highly accurate two-region approximation to N(S) is given by the 
combination of (62) and (77). We have determined, however, that 
(77) alone predicts N(S) with an accuracy of 1 dB or better for S 
between 4.0 and 6.5 for all spectra. Moreover, we have determined that 
N (8S) is at least 119 dB below o? for S = 6.5 and for all spectra. Ob- 
viously, then, (77) can be used over all S for which N(S) is not negli- 
gibly small. We therefore present (77) as our final expression for slope 
overload noise power, with a, --- a, defined by (68) through (71) 
and accurately approximated by (72) through (75) in virtually all 
cases. 


VII. NUMERICAL RESULTS 


The new result has been applied to the two families of spectra shown 
in Fig. 12. The Butterworth spectrum is characterized by an upper 
truncation frequency (W) and the ratio (@1) of the 3-dB corner fre- 
quency to W. Note that, when 8; > ©, the spectrum approaches that 
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TasLe IJI—Np/N, IN DB FoR THE BUTTERWORTH SPECTRUM 





a 0.036 0.111 0.289 0.556 
(61 = 0.02) (61 = 0.068) (61 = 0.25) (6, = ~) 
1.0 +2.85 +1.76 41.44 41.70 
2.0 +0.41 —0.20 —0.22 —0.10 
3.0 +0.43 —0.01 —0.04 +0.22 
4.0 +0.96 +0.74 +0.71 +0.81 
5.0 +1.14 +1.12 +1.11 +1.13 
6.0 +0.90 +0.85 +0.80 +0.70 








of bandlimited white noise. The bandpass uniform spectrum is charac- 
terized by an upper truncation frequency (W) and the ratio (82) of 
the lower truncation frequency to W. Note that 8. = 0 corresponds to 
bandlimited white noise in this case and that, as 8. > 1, the spectrum 
becomes infinitely narrowband. 

For each of these two families, the formula coefficients, (a1, @2, a3, G4), 
have been computed as functions of the respective 6-parameter. The 
results for the Butterworth spectrum are given by Table I and those 
for the uniform spectrum by Table II. The values shown are rounded 
to the decimal accuracy required. 

Curves of N ,/c? versus S for several Butterworth spectra of practical 
interest and for the narrowband spectrum are shown in Fig. 13. (The 
curve for the narrowband case represents an upper bound on N(S)/o? 
for all spectra.) Also, an informative closeup that treats more spectra 
but over a reduced range of S is given by Fig. 14. 

Since the new formula is highly accurate for all spectra, it can be 
used to estimate the accuracy of the Protonotarios formula, (10). 
Table III compares N p(S) with N ,(S) for various Butterworth spectra 
spanning the range of y = bj/babo of practical interest, while Table 
IV gives the comparison for various two-band spectra. The larger 
errors in Table IV must be regarded as untypical, since the two-band 
spectrum itself is rather artificial. For more typical spectra, we can 


TaBLE IV—N p/N, IN DB FOR THE Two-BAND SPECTRUM 





’ 
xX 0.111 0.289 0.556 1 
1.0 +5.04 +3.55 +2.49 +1.52 
2.0 +2.45 +1.46 +0.75 +0.06 
3.0 +1.97 +1.23 +0.70 +0.19 
4.0 +2.05 +154 +1.12 +0.71 
5.0 +1.79 +1.55 +1.30 +1.02 
6.0 +0.84 +0.88 +0.74 +0.50 
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expect the errors in the Protonotarios formula for given S and y to be 
more like the entries in Table III. 


VIII. CONCLUSION 


Aside from providing a new formula for slope overload noise power, 
this work has served to resolve the disparities between previously 
published formulas and can be used to identify the accuracies and 
limitations of those most widely used. Furthermore, the new result 
and the underlying methods of analysis can be extended to treat other 
important topics that have been mostly ignored before. Among these 
are the average duty cycle of slope overload, the slope overload noise 
spectrum, the effects of leaky feedback integrators on slope overload 
noise power, and the slope overload noise power for certain non- 
Gaussian processes. 
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Tube Waveguide for Optical Transmission 


By D. MARCUSE and W. L. MAMMEL 
(Manuscript received September 25, 1972) 


The dielectric optical waveguide described in this paper has an annular 
cross section, the refractive index of which is higher than the indices of 
the material inside and outside of the ring. The solution of the eigenvalue 
problem is an approximation which is valid for small refractive index 
differences of the three media. The resulting approximate eigenvalue 
equation is far stmpler than its exact counterpart. The cutoff conditions 
of the first three modes and the eigenvalue of the lowest-order mode are 
presented graphically. 


I. INTRODUCTION 


Dielectric optical waveguides have become interesting since it has 
been demonstrated that the losses of dielectric materials can be quite 
low.'~ The conventional optical fiber waveguide consists of a solid 
dielectric core surrounded by a dielectric material with lower refractive 
index. This structure has the advantage of being particularly simple. 
In this paper we analyze a different type of dielectric optical waveguide. 
Our structure consists of a dielectric tube which is filled and surrounded 
by dielectric material with lower index of refraction. Such structures 
have been analyzed before.t A complete discussion of and references to 
the literature can be found in Ref. 5. However, the exact analytical 
treatment is extremely complicated so that the resulting theories are 
hard to apply to practical situations. Recently A. W. Snyder® has 
shown that the analysis of the conventional optical fiber waveguide 
becomes much simpler if one assumes that the difference of the refrac- 
tive indices of core and cladding material is only slight. D. Gloge’ 
has developed Snyder’s theoretical work even further. On the basis of 
this simplified method of analysis, it is possible to treat rather com- 
plicated optical waveguide structures more simply and still obtain 
results that are applicable to most cases of practical interest. The 
simplification that results from the assumption of only slight index 
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differences is particularly appropriate since most practical dielectric 
optical waveguides satisfy this requirement. 

An approximate eigenvalue equation for the modes of the tube 
waveguide and an equation describing the cutoff conditions for these 
modes are derived in this paper. The connection of the tube waveguide 
with the dielectric slab waveguide in the limit of large tube radii and 
with the solid-core optical fiber in the limit of small tube radii is 
discussed. The cutoff condition for the first three modes is represented 
graphically and the eigenvalue of the first mode is plotted. 

For simplicity we assume throughout that the index of the material 
inside of the tube is smaller than or equal to the index of the material 
surrounding the tube. 


II. APPROXIMATE MODE ANALYSIS OF THE TUBE WAVEGUIDE 


The geometry of the tube waveguide is shown in Fig. 1. Following 
the ideas of Snyder and Gloge*®’ we express the electric and magnetic 
field components in rectangular cartesian coordinates and write® 
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Fig. 1—Cross section of the tube waveguide. 
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H 1 (6 oH, =) (3) 
oan Ke Ox = Oy 
J 
H t (s 0H, 4 =") (4) 
_— K? oy ~ Ox 
J 


It is assumed that the time and z dependence of the field are determined 
by the factor 
ei(wt—Bz)_ (5) 


8 is the propagation constant in z direction, w is the radian frequency. 
The parameter K; is defined as 


K,=njk? - 8? j =1,2,0r3 (6) 
with 
ke? = wreoto. (7) 
Equations (1) through (4) are exact and are simply four of Maxwell’s 
six equations written in a different form. If the refractive indices are 
all nearly the same, 21 ~ N2 ~ n3 ~ n, we must have 


B = nk, (8) 
so that 

K; <8. (9) 
It is thus apparent that the transverse field components are much 
larger than the longitudinal components if the refractive indices are 
all very nearly the same. The longitudinal field components are 
obtained as solutions of the wave equations.* We use the expressions 

1A;K; 


ie oars {Z,41(K,r) sin (v + 1)¢ + Z,_1(K;r) sin (v — 1)¢} (10) 
n; 


€o tA ;K; 
é \< J (ZreKp) 008 ( + 1)6 
as — Z,(K) cos (v — 1)¢}. (11) 





H, 


The exponential factor (5) has been suppressed. Z,(K,r) is a cylinder 
function and A; is an arbitrary amplitude factor. Substitution of (10) 
and (11) into (1) through (4) results in 


E, = A,Z,(K;r) COs vd, (12) 
A,=-—- nsdn] Z,(K;r) COS vd, (13) 
Ko 
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and 
E,=0, 4H,=0. (14) 


These equations are approximations that hold provided that (8) is 
applicable. 

In the three regions shown in Fig. 1, we use the following cylinder 
functions and parameters, assuming 71 > n2 > Ns, 


2 
—ikKs; = 6 = (6? — nzk?)? 
ae 8 oe for OSr<b (15) 
Z, = J,(i6r) 
2 
Kye = 6h) 
ee es a for <r<a (16) 
Z, = J,(xr) + BN,(«r) 
2 
BH ay ES hk 
— = 6 a 4 for aSrso. (17) 
Z, = H, (tyr) 


J,, N, are the Bessel and Neumann functions and H” is the Hankel 
function of the first kind. We use the functions with imaginary argu- 





0.001 2 4 6 8901 2 4 6 8914 2 4 6 8 | 
d/a 


Fig. 2—The cutoff value of the normalized frequency V for G = 0 as a function 
of the ratio of wall thickness d to tube (outer) radius a. The cutoff value for the 
mode v = 0is V = 0. 
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Fig. 3—The cutoff value of V as a function of d/a for G = 0.7. 


ments instead of modified functions. This practice is in agreement with 
Jahnke and Emde® as well as with Gradshteyn and Ryzhik.” The 
amplitude coefficient A; appearing in (10) through (13) also assumes 
different values in the three regions of space, 7 = 1, 2, 3. 

In order to be able to match the boundary conditions, we transform 
the cartesian components to vector components in cylindrical polar 
coordinates 


EH, = — E,sing + E,cos¢ 
= 24;Z,(K,r)Lcos (v + 1)¢ + cos (v — 1)¢] (18) 
H,= piAn|— Z.(K,r)[sin (v + 1)¢ — sin (v — 1)¢]. (19) 
Ko 


The boundary conditions require that E., Ey, Hz, and Hy are con- 
tinuous at r = b and at r = a. The boundary conditions thus provide 
us with eight equations. However, our approximate treatment of the 
modes yields only four arbitrary constants, Ai, A», A3, and B. The 
situation is further aggravated if we observe that the field components 
contain the functions sin(v + 1)¢ and sin (v — 1)¢ and corresponding 


428 THE BELL SYSTEM TECHNICAL JOURNAL, MARCH 1973 


expressions with the cosine function. Since the boundary conditions 
must hold for all values of ¢, we must equate the coefficients of the 
functions with y + 1 independently of those with » — 1, thus doubling 
the number of equations. In spite of this apparently hopeless situation, 
approximate solutions are possible. The equations that result from 
the requirement that the tangential components of E are continuous at 
the boundaries differ from the corresponding equations resulting from 
the boundary conditions for H only by factors n;. Since our approach 
is based on assuming that all n; are nearly the same, we set these 
factors all equal to the same average value n and, after dividing by n, 
obtain equations that duplicate those obtained from the boundary 
conditions for E. To the approximation considered here, the E and H 
boundary conditions lead to the same set of equations, reducing their 
total number to one half of the original number. To satisfy the con- 
tinuity of H, and Ey at r = b, we equate the coefficients of sin (v + 1)¢ 
and obtain 


_ Jo(xb) + BN, (xb) 


A; 
J,(i6b) 


Ay (20) 


and the determinantal condition 
: J y41(26b) 7 J y4a(xb) + BNy41(xb) 


one (21) 
J (i0b) J(xb) + BN,(xb) 


Equating the coefficients of sin (v — 1)¢, we obtain a similar result 
with the only difference that » + 1 is replaced by v — 1 in (21). 
Equation (20) remains the same. With the help of the recursion 
relations? for cylinder functions, it is easy to prove that eq. (21) is 
not changed if »y + 1 is replaced by v — 1. We thus find ourselves in 
the happy position of being able to satisfy—at least approximately—the 
eight boundary conditions at r = b with the two available coefficients. 
In an analogous fashion, we obtain from the boundary conditions at 
r = a the equation 
J (xa) + BN,(xa) 


A, = A 22 
eae (22) 





and the determinantal condition 


Hoy x(tya) Sana) + BN 41(Ka) 


ue H®(iya) ™ UEeay BN Cea) v2) 
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0.01 





Fig. 4—Same as Fig. 3, G = 1.5. 


The eigenvalue equation is obtained by eliminating B from (21) and 
(23). We find 
KJ (40d) J 4.1(kb) — 107 ,(xb) J ,41(20b) 
KJ (100) Nyyi(kb) — t0N,(xb)J,41(20b) 
KH, (¢ya)Jy41(ka) — tyJ,(xa)H,41(éva) 
© eH Eya)Noaa(ca) — FN, (ea)H (ya) 








(24) 


The constant B is given by 


Bik. ee) _ cls (25) 
kJ ,(70b) Ny13(xb) — 16N,(xb)J,41(20b) 





It can be shown that the eigenvalue equation (24) specializes to the 
eigenvalue equation (d = a — b) 


= K(y + 8) 
K2 — ¥@ 


tan xd (26) 


of the asymmetric slab waveguide in the limit a— © and b >. The 
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distinction between TE and TM modes of the slab waveguide is lost 
in our approximation. 

We are mainly interested in the cutoff conditions of the modes in 
order to find the separation between the lowest and the next higher 
mode which determines the range of single-mode operation. Cutoff is 
defined by the condition 


y= 0. (27) 


In order to determine the cutoff frequency, we use (27) in (24) and 
obtain the cutoff equation 


b b b b 
Jy (IGV —) Ivar ( V—) — iGd, | V — ) Joni ( 1GV - 
d d d d 
b b b "24-2 
J iGV — Ny4i V- — iGN, V- Ja iGV — 
d d d d 


= 28) 


The parameters appearing in this equation are defined as follows: 





V = (ni — n3)'kd (29) 
and 
2 a2, 
Ne — N3\? 
G = (“ = ~) . (30) 
1 2 


For future reference, we derive an interesting relation for the mode 
with v = 0. The cutoff condition (28) allows the solution V = 0 for 
v = 0. We find as the condition that V = 0, the relation 


2 Dg 

d Ni — No\? 1 

-=1- = 1—-—_—__—.. (31) 
(1 + G2)! 





In the limit of infinite radii, a— «© and b—, we obtain from (26) 
the cutoff condition for n; > nz = nz 


V = arctan G. (82) 
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0.001 
Fig. 5—Same as Fig. 3, G = 5.3. 


III. DISCUSSION AND NUMERICAL RESULTS 


We restrict our discussion to the case nz 2 n3. The cutoff values of 
V [see eq. (29) ] as functions of d/a are plotted in Figs. 2 through 6 
as solid lines for the first three modes, v = 0, »y = 1, v = 2. The V 
parameter characterizes the tube waveguide in the limit of large radii 
where it can be regarded as an asymmetric slab waveguide that is 
bent into a circle. For large values of d/a, the waveguide approaches 
an optical fiber with solid cylindrical core. Such a fiber is characterized 
by a different V parameter which can be expressed in our notation as 
(a/d)V. In the limit d/a = 1, both parameters coincide. The parameter 
(a/d)V is plotted in the figures as a dotted line only for the mode 
vy = 1]. Each curve, Figs. 2 through 6, is plotted for a different value 
of G [see eq. (30) ]. The relation between the G values and the ratios 
of n1/nz and n1/n3 is shown in Fig. 7. 

Since the tube waveguide approaches the asymmetric slab in the 
limit d/a = 0 and the solid-core fiber in the limit d/a = 1, the cutoff 
values V = V. can be predicted at these limits. For d/a = 1, the 
cutoff value of the solid-core fiber is obtained from 


Jy31(Ve) = 0. (33) 
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The values of V, at d/a = 1are, therefore, V. = Oforv = 0, V. = 2.405 
forv = land V, = 3.83 fory = 2. 

At d/a = 0, we find the cutoff values for V from the asymmetric 
slab formula (32). 

We see from Figs. 3 through 6 that the curve with v = 0 seems to 
end on the horizontal axis. Actually, the curve must be continued along 
the horizontal axis at V = 0 to the point d/a = 1. However, to the 
left of the point where the v = 0 curve reaches the horizontal axis, 
V = 0 is not a legitimate solution of the cutoff equation. The d/a 
value at which the v = 0 curve touches the horizontal axis is given by 
(31). This phenomenon finds the following simple explanation (private 
communication from E. A. J. Marcatili). As the mode approaches its 
cutoff value, we have y = 0, and also @ = 0, and V = 0. There is no 
longer any transverse field variation so that the field ‘‘senses’’ only 
the average value 7? of the dielectric constant of the tube, 


xb?n3 + x(a? — b2)n; 


Tra? 


72 = (34) 


The mode is no longer guided if the average dielectric constant 7? is 





2 4 6 8 0.01 2 4 6 8 0.1 2 4 6 8 1 
d/a 


Fig. 6—Same as Fig. 3, G = 9.6. 
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Fig. 7—The parameter G is shown as a function of (m1/n3 — 1) for several values 
of n/N. 


equal to the dielectric constants nj of the outer medium. The condition 
n? = ni also leads to (31). 

Figures 2 through 6 show that the cutoff points for all three modes 
(vy = 0, 1, and 2) approach each other arbitrarily closely for small 
values of d/a. It is thus apparent that single-mode operation in the 
mode v = 0 becomes increasingly more difficult to achieve as the 
ratio d/a is decreased. 

Finally, Fig. 8 shows a plot of the transverse decay parameter ya 
as a function of d/a. The curve is intended to provide information 
about single-mode operation of the tube waveguide. However, since 
it is desirable to operate with a tightly guided mode field, we chose as 
the operating point for V the cutoff value of the second mode, V = Vai. 
It is a remarkable feature of the curves of Fig. 8 that they become 
independent of d/a for small values of this ratio. However, from the 
point of view of field confinement, we must remember that it is not 
ya but rather yd that characterizes the rate at which the field decays 
from the guiding structure—the tube—in radial direction. The crosses on 
the curves in Fig. 8 indicate the points at which we have yd = 1, the 
circles indicate the points yd = 0.5. 
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Fig. 8—This figure shows the parameter ya as a function of d/a for several values 
of G. V is taken equal to the cutoff value Va of the mode v = 1. 


IV. CONCLUSIONS 


The eigenvalue equation for the modes of the tube waveguide has 
been derived with the help of an approximate technique. The cutoff 
values of the normalized frequency parameter V [defined by (29) ] are 
presented for the first three modes of the tube waveguide and the 
transverse decay parameter ya is plotted for the lowest-order mode 
under the assumption that this mode is operated at the point of cutoff 
of the next-higher-order mode. 

One may wonder if it is possible to use the tube waveguide with 
single-mode operation to increase the mode radius compared to the 
mode radius of the lowest-order HE,; mode of the conventional fiber. 
This question has been studied under the assumption that both wave- 
guides provide an equal amount of field confinement. The criterion for 
field confinement is the radial decay parameter y. It was consequently 
assumed that the y values of both waveguides are identical. Surpris- 
ingly, it was found that under the condition of single-mode operation 
for each guide and the requirement of equal field confinement, the 
mode radius of the tube waveguide is only slightly (approximately 
40 percent) larger than the radius of the HE, mode of the optical 
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fiber. A large mode radius appears desirable from the point of view 
of alignment tolerances of fiber splices. With a large mode radius, the 
alignment of the waveguide centers would be less critical. The tube 
waveguide does not appear to offer a significant advantage for this 
purpose. 
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The Interrupted Poisson Process 
As An Overflow Process 


By ANATOL KUCZURA 
(Manuscript received September 18, 1972) 


Traffic overflowing a first-choice trunk group can be approximated 
accurately by a simple renewal process called an interrupted Poisson 
process-a Poisson process which is alternately turned on for an ex- 
ponentially distributed time and then turned off for another (independent) 
exponentially distributed time. The approximation is obtained by matching 
erther the first two or three moments of an interrupted Poisson process to 
those of an overflow process. Numerical investigation of errors in the 
approximation and subsequent experience has shown that this method of 
generating overflow traffic 1s accurate and very useful in both simulations 
and analyses of traffic systems. 


I. INTRODUCTION 


A Poisson process which is alternately turned on for an exponentially 
distributed time and then turned off for another (independent) 
exponentially distributed time will be called an interrupted Poisson 
process—it can be viewed as a Poisson process modulated by a random 
switch. It was suggested by W. S. Hayward of Bell Laboratories that 
such a process be used to simulate overflow traffic. We will show that 
the interrupted Poisson process provides a simple and accurate method 
of simulating overflow traffic. 

The objective is to reduce the cost of computer simulations of 
traffic systems by using the interrupted Poisson process to model the 
overflow traffic. Generation of actual overflow traffic by simulating 
the behavior of the trunk group from which it overflows is time con- 
suming since a record must be kept of all calls which are offered to the 
subtending trunk group, whether they contribute to the overflow 
traffic or not. This is especially true when the traffic is overflowing a 
large trunk group. Moreover, the interrupted Poisson process provides 
a simple, approximate description of the overflow traffic and conse- 
quently facilitates analytical studies. 
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This method of generating overflow traffic has been used successfully 
in many studies of traffic systems. Examples of its application can be 
found in Refs. 1, 2, and 3 and also in numerous unpublished works. 
Since the method has been found to be very useful and requests for 
wider dissemination have been received by the author, this paper has 
been prepared. 

In Section II we derive the distribution of the number of busy 
trunks in an infinite trunk group when the offered traffic is generated 
by an interrupted Poisson process. The corresponding distribution for 
an overflow input has been computed by Kosten.4 Now, matching the 
first three moments of the two distributions, we obtain equations for 
the parameters of the interrupted Poisson process. We also give these 
equations for a two-moment match. The errors committed in approxi- 
mating the distribution given by Kosten are also examined. In Section 
III we derive the interarrival time distribution for a traffic stream 
generated by an interrupted Poisson process. 


II. MOMENT-MATCH EQUATIONS 


Let the interrupted Poisson traffic be offered to an infinite trunk 
group. Let \ be the intensity of the Poisson process, 1/7 be the mean 
on-time of the random switch, 1/w be the mean off-time, and 1/y be 
the mean service time. Let the state of the system be described by 
(m, n) with state probabilities p(m, n) where m is the number of 
servers busy, and n is the state of the switch taking on the value of 1 
or 0 according to whether the process is on or off. 

The equilibrium equations for the stationary state probabilities are 


(mp + «)p(m, 0) = yp(mn, 1) + (m+ Lup(m + 1,0), m=O, 
(mp + y + A)p(m, 1) = wp(m, 0) + (m + l)up(m + I, 1) 
—  +hp(m~ 1,1), m2l, (1) 
(y + )p (0, 1) = wp(0, 0) + up (1, 1). 
To solve this system of equations we introduce the probability generat- 
ing function 


G(z) = L p(m)z = Giz) + Go(2), 
where . 
Gi(z) = L p(m, Lz", Go(z) = L pm, 0)”, 
and . 
p(m) = p(m, 1) + p(m, 0). 
Note that G:(1) = >} p(m, 1) is simply the probability of the switch 
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being on and is given by w/(y + w). Similarly, Go(1) = > p(m, 0) 
= y/(y +) is the probability of the switch being off. We will now 
derive the differential equations for G, and Gp» and obtain their 
solutions. 

From (1) and the definition of G(z) we have 


u(z — 1)Go(z) + wGo(z) — yGi(z) = 0, 
uz — 1)G4(z) + (vy +A — Az)Gi(z) — wGo(z) = 0. 
This system of equations is coupled. At the price of increasing the 


order, we can decouple the system by means of differentiation and 
simple substitution: 


ulz — 1)G0 (2) + Te ty tw — Az — 1) G02) 
— ® Gale) =(0, (2) 
m 
u(z — 1)GY (2) tLe +y +o — Az — 1) IG) 
nN 
a= (ab wea = 6. 
m 


Since we are interested in the moments of the distribution of the 
number of busy servers, it will be convenient to obtain solutions of (2) 
valid about z = 1. Subsequently, the series can be rearranged at the 
origin to yield the state probabilities. 

The change of variable 


rd 
eS =e 1) 
be 


transforms (2) into 


EQo(£) + (e — E)Qo(E) — BQo(E) = O, 





(3) 
EQ1(é) + (e — &)Qi(E) — (1 + B)Qi(E) = 0, 
where 
e= 1 oP u aE 2 ? B = ay 
ue Ki 
and 


Q(t) = Gi), Q2(f) = Ge(z). 
Equations (3) have the solutions 
Qo(é) = CiF1(B; €; &), 
Qi(é) = Dif i(l + B; ; &), 
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where 
RG 
a;b;c) = 
ce n=0 (b) nn! 
with 


(q)n =a(at+ I@+2)---@tn—-1), @o=1, 


is the confluent hypergeometric function and C and D are arbitrary 
constants. The conditions G; (1) = w/(y + w) and Go(1) = y/(y + w) 
determine C' and D and hence 


oa [a a 1)| 
z) = W1| Bye;-(2 — 
Y +o LU 


+ 
y+o 


The power series in (2 — 1) on the right-hand side of (4) is con- 
vergent since ,f; is an entire function. The factorial moments of the 
distribution of the number of busy servers, being the coefficients of 
the expansion (4), are simply 

(w)n 
Gel) = n = 0, 1, 2, oe Ey (5) 
(y 5S w) n 

where we have made the normalization » = 1. For computational 

purposes, the following recurrence relation is useful : 

wtn 
G@+n(1) = OE) cee (6) 
(y+@+n) 

To obtain the state probabilities p(m) = p(m,0) + p(m,1) we 
rearrange the series (4) at the origin and identify the coefficients in 
this new expansion as the state probabilities 








Fy E +856-@ = 1). (4) 


(= a 
p(m) = — aoe GO) 
m! k=0 
d)a-™ @) 5 
I RE) M4 (7) 


=p (j—mly ae 


We now show that the interrupted Poisson process provides an 
accurate method of generating overflow traffic. Let a erlangs of Poisson 
traffic be offered to an Erlang B system of c trunks. Let the overflow 
traffic be routed to an infinite trunk group and let Y be the number of 
busy servers in the infinite trunk group under statistical equilib- 
rium. The factorial moments 1/(,) of Y and the state probabilities 
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f(m) = pLY = m] have been computed by Kosten :4 








Wes Zia ahe) 
on(€) 
(8) 
qetm ce) (—a)* 
f(m) = —— 
clm! k=0 Kon 4m(C) 
where 
ac fj i— iy at 
g(c) =—, oi(e) = = ( 7 \— j=1,2,-°-. 
c! = 1 (ec — 7) 


For a more accessible reference which gives the derivation of the 
factorial moments, see the appendix prepared by J. Riordan in Ref. 5. 

Now consider an interrupted Poisson traffic of original intensity \ 
offered to an infinite trunk group, and let X be the number of busy 
servers under statistical equilibrium. The factorial moments of X and 
the state probabilities are given by (5) and (7) respectively. In this 
system we have three parameters, A, y, and w, which are to be chosen 
so that the interrupted Poisson process gives the best approximation 
to the overflow process. In the present analysis we choose the moment 
approximation and, in particular, take the factorial moments. Thus we 
require that 


God) = Ma), n= 1, 2,3, (9) 

and define the error 
E, (a,c) = f(n) — p(n), n=0,1,2,---. (10) 
With the aid of the recurrence relation (6), we can express (9) as 
(——) aie. = 00102) (11) 


where 
M (n+1) on(C) 


aM (ny @ngalc) 
Equations (11) have the solution 
5o(61 — 50) — So(d2 — 61) 
(61 — 80) — (62 — 41) 


60 A- a6, 
ga (—"), (12) 
rN \d1 — 50 


“( _ ~*) 
yur : 
a 60 


) 
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The parameters \, w, and y can be computed from (12) whenever 
a and c are specified. Often, however, one is concerned with final-route 
traffic in which only the mean, a, and variance, v, (or peakedness ratio 
z = v/a) are known. There are two ways to proceed in this case. 

First, using Wilkinson’s equivalent random method,® determine S, 
the number of trunks, and A, the equivalent random load correspond- 
ing to the overflow traffic of mean a and variance v. Now set a = A 
and c = S and use eqs. (12) to compute X, w, and y for a three-moment 
match. 

A second way to proceed is to determine the equivalent random 
load A as before, and set \ = A in the last two equations of (12) to 
compute w and y for a two-moment match. A satisfactory value of the 
equivalent random load A is given by Rapp’s approximation :* 


A = az + 82(z — 1), (13) 


where z is the peakedness ratio v/a. In terms of a, z, and A, these 
equations can be written as 


a {/A-—-a 
w= =( -1), 
A\z-1 


(Ep 


Note that for a two-moment match it is necessary to fix one of the 
parameters A, w, or y. However, for a positive solution, they cannot 
be chosen arbitrarily—for positivity we must choose \ > a, such as 
in (18). 

To illustrate the procedure, let us take an example with overflow 
traffic of mean 0.61 and variance 0.95. By Wilkinson’s equivalent 
random method, we obtain S = 4 trunks and A = 8 erlangs. The 
first three moments from (8) with c= S = 4 and a= A = 8 are 





(14) 


May = 0.618 
M2 = 0.708 (15) 
M 33) = 1.025. 


The three-moment match yields an interrupted Poisson process with 
parameters 


A = 2.553 
w = 0.646 
y = 2.022, 


and of course the same first three moments as in (15). 
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For the two-moment match, we set \ = A = 3 and from (14) obtain 
w = 0.669 
y = 2.621. 


The first two moments will be equal to Ma) and M) of (15) and the 
third moment is found to be 


G®(1) = 1.049 


which is not significantly different from M3) in (15). Computing the 
state probabilities, we obtain the following result in which the two- 
moment match of the negative binomial fit® was included for 
comparison: 


Exact State Three-Moment Two-Moment Negative 


State Probabilities Match Match Binomial 
0 0.6164 0.6161 0.6149 0.6086 
1 0.2312 0.2318 0.2344 0.2464 
2 0.0965 0.0962 0.0953 0.0924 
3 0.0372 0.0372 0.0366 0.0337 
4 0.0130 0.0130 0.0129 0.0122 
5 0.0041 0.0041 0.0042 0.0043 
6 0.0012 0.0012 0.0012 0.0015 
7 0.0003 0.0003 0.0003 0.0005 
8 0.0001 0.0001 0.0001 0.0002 


We now examine the error FE, (a,c). Since little additional com- 
putation is required to obtain the three-moment match, we feel it 
should be done to obtain a better fit. Consequently, we examined the 
error for a three-moment match only. Calculations of f(n), the state 
probabilities as they are given exactly, and p(n), the interrupted 
Poisson approximation, have been made. This was done for groups of 
1, 2, 4, 8, 16, 32, 40, and 48 trunks in the primary group with offered 
occupancies of 0.75 and 1, and for a group of 64 with offered occupancy 
of 0.75. For larger trunk groups, significant loss of accuracy prevented 
successful computation. We note that for the case of one trunk the 
approximation is exact. Where comparison could be made with pre- 
viously reported results for the negative binomial and the confluent 
hypergeometric approximations,’ it was found that the three-moment 
match using the interrupted Poisson process gave uniformly better fit 
to the state probabilities. 

A typical result of the computations made is displayed in Table I. 
It was found that for a fixed occupancy the errors would increase, 
reach a maximum, and then decrease as the number of trunks was 
increased. This behavior can be seen in Tig. 1 where 


E = max |£,(a,c)| = max | f(n) — p(n)| 
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TasLeE J—TsHree—-Moment Matcu 
(c = 16, a/c = 0.75, \ = 6.83, w = 0.366, y = 3.09) 











E, (12,16) E, (12,16) 
State n f(n) p(n) = f(n) — p(n) f(n) 
0 0.64866359 0.64765765 0.00100593 0.00155078 
1 0.17173954 0.17395867 —0.00221912 —0.01292144 
2 0.08460030 0.08381668 0.00078362 0.00926265 
3 0.04523738 0.04463748 0.00059990 0.01326127 
4 0.02430285 0.02417444 0.00012841 0.00528357 
5 0.01280780 0.01290307 —0.00009527 —0.00743819 
6 0.00655960 0.00668314 —0.00012354 —0.01883330 
7 0.00325188 0.00333032 —0.00007845 —0.02412386 
8 0.00155792 0.00158885 —0.00003093 —0.01985106 
9 0.00072097 0.00072378 —0.00000280 —0.00388859 
10 0.00032234 0.00031441 0.00000793 0.02460252 
1] 0.00013930 0.00013020 0.00000910 0.06530564 
12 0.00005822 0.00005142 0.00000681 0.11694969 
13 0.00002356 0.00001937 0.00000418 0.17756000 
14 0.00000923 0.00000697 0.00000226 0.24472624 
15 0.00000351 0.00000240 0.00000111 0.31590026 
16 0.00000129 0.00000079 0.00000050 0.38876238 
17 0.00000046 0.00000025 0.00000021 0.46060149 
18 0.00000016 0.00000008 0.00000009 0.53013187 
19 0.00000005 0.00000002 0.00000003 0.59587179 














is plotted against c for different values of the offered occupancy a/c 
and again in Fig. 2 where the relative error corresponding to the state 
n at which |£,,(a, c)| attained its maximum is plotted. 


II]. INTERARRIVAL TIME DISTRIBUTION 


In analytical studies of systems, a description of overflow traffic is 
sometimes needed in terms of the interarrival time distribution.? This 
distribution is complicated and may be difficult to compute whenever 
the size of the trunk group which the Poisson traffic is overflowing is 
large.® If the interrupted Poisson traffic is used to generate the overflow 
traffic, the resulting interarrival time distribution, say A (¢), is simple. 
Indeed, we show that this distribution is given simply by the mixture 
of two exponential distributions : 


A(t) =k(l — eo) + k(1 — e-*), (16) 
m= 7A totyt VA+et 7)? — dro}, 
r= 7A toty- VAtowt 7)? — Aro}, 


where 


\— 12 
ky = ) 
T1 — Te 


kp = 1 — ky. 
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Fig. 1—Maximum absolute error E. 


The proof goes as follows. Let W, be the waiting time from ¢t = 0 
until the time of the nth arrival and H,(¢) be the distribution of W,. 
To obtain the interarrival time distribution, it is not necessary to 
find H, (t) for all n. The distribution H,(¢) and proper choice of initial 
conditions at ¢ = 0 is sufficient to find A(t). We include the more 
general case here for completeness. 

If N() counts the number of arrivals in (0, t) and 


px(t) = PLN (é) = k, 
then 
H,(t) =1— 2 pilt). (17) 
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Fig. 2—Relative absolute error corresponding to Fig. 1. 


This equation follows from the observation that W, > t if and only 
if N@) <n. 

Taking the Laplace-Stieltjes transform of both sides of (17) we 
obtain 


an(s)=1—s > wx(S), (18) 
where oT 


an(S) 


I 


i e “dH, (t), 
0 


w(s) = [ e*'p,(t)dt. 
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The initial conditions used to obtain (18) are 
Ts k=0 


p(Ot) = 0, pe. (19) 


We will compute 7; (s) and hence determine a, (s). From the expression 
for a»(s) we then determine the interarrival distribution A (f). 

Let pim(t) be the probability that there were & arrivals in (0, 0), 
given that an arrival occurred at ¢ = 0 and that at the instant ¢ the 
switch is on if m = 1 and off if m = 0. These functions satisfy the 
system of differential equations 


por(t) = wpoo(t) — (A + y)poi(d), 
pratt) = wpro(t) = (A ae ¥)per(t) + Apr-1,1(t), k= i; 2, adi) (20) 
pro(t) > wp ro(t) + ypx(t), k= 0, 1; 2, Tins 2) 


and the initial condition poi (0) = 1. 
Taking the Laplace transform of (20), we obtain 


S7o1(s) = w7oo(s) = (v -+- )701(8) + 1, 
Smxi(s) = waeo(s) — (A + y)me1(s) + Ame-1,1(8), k= 1,2, -+<, 
st x0(s) = — wrxo(s) + Ar E1(S), k= 0, I, 2, Bnet 


where 
Wxj(S) -{ e*'nxj(t)dt. 
0 


This system of difference equations can be solved for rz0(s) and 
m11(8) and hence 7; (s), since 7;(s) is the sum 7z0(s) + 721(s). Omitting 
the details, we have 
eer) 


——— =0,1,2,---, (21 
7(s) 7(s) oe 


11(S) => 

where 
Tis) = PF +A+t+yt+to)s + drow. 

Substituting (21) into (18), we obtain 


A(s + w)7]” 
no roe 


The interarrival time distribution is given by the distribution of the 
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waiting time until the first arrival and hence its Laplace-Stieltjes 
transform is given by a:(s). Inverting a:(s), we obtain eq. (16). 
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